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Abstract 


The  work  presented  in  this  document  contributes  to  the  ROBR  (Reconfigurable  Omni  Band 
Radio)  project  started  by  the  Defence  Research  Establishment  Ottawa  and  the 
Communication  Research  Centre  in  1997.  ROBR  is  a  testbed  implementation  of  a 
reconfigurable  satellite  communications  (satcom)  terminal  that  meikes  xise  of  a  software 
communications  architecture.  Such  a  system  can  enable  the  use  of  a  single  groimd  terminal  to 
communicate  over  multiple  satellite  communications  or  terrestrial  links  by  supporting 
multiple  standards.  The  ROBR  hardware  architecture  includes  a  microprocessor  and  several 
digital  signal  processor  (DSP)  boards.  The  objective  of  this  report  is  to  document  the  work 
done  to  provide  a  set  of  reconfigurable  digital  filters  for  use  in  the  ROBR.  Five  infinite 
impulse  response  (HR)  filtering  modules  and  four  finite  impulse  response  (FIR)  filtering 
modules  have  been  implemented.  The  fimction  of  these  modules  is  to  compute  the 
coefficients  of  a  desired  filter  design.  Also,  HR  and  FIR  signal  processing  modules  have  been 
implemented  to  process  digital  signals  using  the  computed  coefficients.  The  modules  have 
been  implemented  in  the  C  programming  language  and  are  targeted  for  use  on  a  DSP  chip. 
The  implementation  of  the  modules  has  been  verified  and  compared  with  the  results  obtained 
with  the  Signal  Processing  toolbox  from  MATLAB. 

R^SUM^ 

Les  travaux  presentes  dans  le  present  document  contribuent  au  projet  ROBR  (radio  omni- 
bande  reconfigurable)  entrepris  par  le  Centre  de  recherches  pour  la  defense,  Ottawa  et  le 
Centre  de  recherches  sur  les  communications  en  1997.  II  s’agit  de  la  mise  en  oeuvre  d’un 
prototype  de  terminal  reconfigurable  de  teleconununications  par  satellite  qui  a  recours  a  ime 
architecture  logicielle  pour  les  communications.  Ce  terminal  pent  permettre  Tutilisation  d’lm 
seul  terminal  au  sol  pour  assurer  des  communications  au  moyen  de  plusieurs  liaisons  de 
communications  par  satellite  ou  de  plusieurs  liaisons  de  Terre  en  vertu  de  plusieurs  normes. 
L’architecture  matMelle  du  ROBR  comprend  un  microprocesseur  et  plusieurs  cartes  de 
traitement  numerique  des  signaux  (DSP).  Le  present  rapport  a  pour  but  de  documenter  le 
travail  effectue  pour  la  foumiture  d’un  jeu  de  filtres  numeriques  reconfigurables  en  vue  de 
son  utilisation  dans  le  ROBR.  Cinq  modules  de  filtrage  des  reponses  impulsionnelles  infinies 
(RII)  et  quatre  modules  de  filtrage  des  reponses  impulsionnelles  finies  (RIF)  ont  ete  mis  au 
point.  Ces  modules  servent  au  calcul  des  coefficients  de  la  conception  ddsiree  des  filtres.  Des 
modules  de  traitement  des  signaux  RII  et  RIF  ont  aussi  ete  mis  au  point  pour  le  traitement 
numerique  des  signaux  au  moyen  des  coefficients  calcules.  Ces  modules,  configures  en 
langage  de  programmation  C,  doivent  etre  utilises  sur  des  puces  DSP.  Le  fonctionnement  des 
modules  a  ete  verifie  et  compare  aux  resultats  obtenus  au  moyen  du  produit  Signal 
Processing  Toolbox  de  MATLAB. 
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Executive  Summary 


The  work  presented  in  this  document  contributes  to  the  ROBR  (Reconfigurable  Omni  Band 
Radio)  project  started  by  the  Defence  Research  Establishment  Ottawa'  and  the 
Communication  Research  Centre  in  1997.  ROBR  is  a  testbed  implementation  of  a 
reconfigurable  satellite  communications  (satcom)  terminal  that  can  support  multiple 
standards.  Such  a  system  potentially  allows  the  Canadian  Forces  to  use  a  single  ground 
terminal  to  communicate  over  multiple  satellite  communications  or  terrestrial  links. 

The  ROBR  hardware  architecture  includes  a  microprocessor  and  several  digital  signal 
processor  (DSP)  boards.  The  objective  of  this  report  is  to  dociunent  the  work  done  to  provide 
a  set  of  reconfigurable  digital  filters  for  use  in  the  ROBR.  Five  infinite  impulse  response 
(HR)  filtering  modules  and  fo\ir  finite  impulse  response  (FIR)  filtering  modules  have  been 
implemented.  These  modules  will  be  integrated  in  the  DSPs  and  processed  by  them.  The 
fimction  of  these  modules  is  to  compute  the  coefficients  of  a  desired  filter  design.  Also,  HR 
and  FIR  signal  processing  modules  have  been  implemented  to  process  test  input  signals  using 
the  filter  coefficients  generated. 

The  implemented  techniques  to  design  HR  filters  are  based  on  transformation  of  continuous¬ 
time  HR  systems  into  discrete-time  HR  systems.  Two  conversion  methods  have  been 
implemented  to  produce  a  discrete  filter  design  from  an  analog  filter  design:  the  “Bilinear 
Transformation”  and  the  “Impulse  Invariance”  methods.  Five  major  types  of  HR  filters  have 
been  implemented,  the  Butterworth  filter,  the  Chebyshev  filter,  the  inverse  Chebyshev  filter, 
the  Elliptical  filter  and  the  Bessel  filter.  The  computation  of  the  coefficients  for  the  HR 
filtering  modules  is  separated  into  three  steps:  the  analog  design  computation,  the  conversion 
of  the  analog  design  into  a  discrete  design  and  the  factorization  step  which  produces  the 
coefficients  of  the  digital  filter. 

The  implemented  techniques  for  computing  the  coefficients  of  the  FIR  filtering  modules  are 
the  “Frequency  Sampling  Design”  method,  the  “Design  by  Windowing”  method  and  the 
“Parks-McClellan”  method  (also  called  the  “Remez  Exchange  Algorithm”).  The  general 
procedure  for  FIR  filter  design  is  to  sample  the  frequency  response  of  a  filter  and  then 
compute  the  inverse  discrete  Fourier  transform  (IDFT). 

An  FIR  implementation  to  compute  the  coefficients  of  a  Gaussian  filter  has  been 
implemented.  The  method  used  to  calculate  the  coefficients  of  this  filter  is  the  “Frequency 
Sampling  Design”  method.  Finally,  a  digital  integrator  was  implemented.  The  function  of  a 
digital  integrator  is  to  sum  a  digital  input  sequence  over  time.  The  Gaussian  filter  module  and 
the  digital  integrator  module  may  be  used  in  the  premodulation  stage  of  a  Gaussian  minimum 
shift  keying  (GMSK)  digital  modulator  for  the  ROBR. 

The  modules  have  been  implemented  in  the  C  language  for  further  implementation  on  a  DSP 
chip.  The  implementation  of  both  HR  and  FIR  filter  modules  have  been  verified  and 
compared  with  the  results  obtained  with  the  Signal  Processing  toolbox  fi-om  MATLAB. 
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The  implementation  of  two  filtering  modules  has  been  successfully  done  on  a  TMS320c6201 
digital  signal  processor  from  Texas  Instruments  mounted  on  a  Daytona  DSP  board  from 
Spectrum. 


Gosselin,  B.,  and  Wilcox,  C.,  2001,  Reconfigurable  Digital  HR  and  FIR  Filters,  DREO  TR 
2001-099,  Defence  Research  Establishment  Ottawa. 
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L’architecture  materielle  du  ROBR  comprend  un  microprocesseur  et  plusieurs  cartes  de 
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travail  effectue  pour  la  foumiture  d’lm  jeu  de  filtres  nummques  reconfigurables  en  vue  de 
son  utilisation  dans  le  ROBR.  Cinq  modules  de  filtrage  des  reponses  impulsionnelles  infinies 
(RII)  et  quatre  modules  de  filtrage  des  reponses  impulsionnelles  finies  (RIF)  ont  6td  mis  au 
point.  Ces  modules  seront  integres  aux  DSP  et  traites  par  eux.  Ils  servent  au  calcul  des 
coefficients  de  la  conception  desiree  des  filtres.  Des  modules  de  traitement  des  signaux  RII  et 
RIF  ont  aussi  ete  mis  au  point  pour  le  traitement  des  signaux  d’entrde  d’essai  au  moyen  des 
coefficients  generes  pour  les  filtres. 

Les  techniques  utilis6es  pour  la  conception  des  filtres  RII  sont  fondees  sur  la  transformation 
des  systemes  RII  a  temps  continu  en  systemes  RII  a  temps  discret.  Deux  methodes  de 
conversion  ont  ete  mises  en  oeuvre  pour  la  production  d’une  conception  de  filtre  discret  a 
partir  de  la  conception  d’rm  filtre  analogique  :  la  methode  de  “  transformation  bilindaire  ”  et 
la  methode  “par  invariance  impulsionnelle  ”.  Cinq  grands  types  de  filtres  RII  ont  ete 
selectionnes  ;  le  filtre  de  Butterworth,  le  filtre  de  Chebyshev,  le  filtre  de  Chebyshev  inverse, 
le  filtre  elliptique  et  le  filtre  de  Bessel.  Le  calcul  des  coefficients  pour  les  modules  de  filtrage 
RII  se  divise  en  trois  etapes :  le  calcul  de  la  conception  analogique,  la  conversion  de  la 
conception  analogique  en  conception  discrete  et  la  factorisation  (qui  produira  les  coefficients 
du  filtre  numerique). 

Les  techniques  selectionnees  pour  le  calcul  des  coefficients  des  modules  de  filtrage  RIF  sont 
la  methode  de  “  conception  par  echantilloimage  de  fi'equences  ”,  la  methode  de  “  conception 
par  fenStrage  ”  et  la  methode  de  “  I’algorithme  d’ optimisation  de  Parks-McClellan  ”  (aussi 
appelee  la  “fonction  de  Remez”).  La  procedure  generale  de  conception  des  filtres  RIF 
consiste  a  echantillonner  la  reponse  en  fi:equence  d’un  filtre,  puis  a  calculer  la  transformee  de 
Fourier  discrete  inverse. 

Une  application  aux  RIF  pour  le  calcul  des  coefficients  d’un  filtre  gaussien  a  ete  effectuee. 
La  methode  utilis^e  pour  le  calcul  des  coefficients  de  ce  filtre  est  la  mdthode  de  “  conception 
par  echantilloimage  de  fi’equences  ”.  Enfin,  im  module  d’integrateur  numMque  a  ete  elabore. 
Un  integrateur  numdrique  sert  au  calcul  de  la  somme  d’une  sequence  d’entree  numerique 
pour  une  periode  donnde.  Le  module  de  filtre  gaussien  et  le  module  d’integrateur  numerique 
peuvent  etre  utilises  a  I’etape  de  la  premodulation  d’lm  modulateur  numerique  de  modulation 
par  deplacement  minimal  avec  filtrage  gaussien  (MDMG). 


Vll 


Les  modules,  configures  en  langage  de  programmation  C,  doivent  etre  utilises  sur  des  puces 
DSP  a  ime  date  ulterieure.  Trois  modules  de  filtrage  ont  ete  implantes  avec  succds  sur  vm 
processeur  numerique  de  signaux  TMS320c6201  de  Texas  Instruments  monte  sur  une  carte 
DSP  Daytona  de  Spectrum.  Le  fonctionnement  des  modules  a  ete  v^rifie  et  compart  aux 
resultats  obtenus  au  moyen  du  produit  Signal  Processing  Toolbox  de  MATLAB. 
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1  Introduction 


With  the  rapid  increase  of  communications  services  and  systems  being  developed  and 
implemented  using  different  standards,  much  effort  has  been  directed  over  the  past  decade  to 
ded  with  issues  relating  to  interoperability  and  compatibility  of  these  various  systems. 
While  global  standardization  can  be  one  solution  to  the  problem,  a  more  practical  remedy 
may  be  to  develop  transceivers  that  can  support  different  frequency  bands  and  different 
waveform  standards  using  a  common  hardware  platform  or  device.  A  feasible  way  to 
achieve  this  is  to  implement  the  waveform  standards  in  software  to  provide  the  required 
flexibility. 

Although  initially  focused  on  the  personal  communications  services  industry,  the  software 
radio  concept  can  be  extended  to  other  applications  such  as  satellite  communications 
(satcom).  Just  as  multiple  standards  are  being  considered  for  integration  on  a  single  device 
(e.g.  cellular  phone),  m^tiple  satcom  waveforms  can  be  implemented  in  software  to  run  on  a 
single  hardware  platform  or  terminal.  As  with  personal  communications  services,  this 
approach  potentitdly  offers  benefits  for  interoperability,  ease  of  future  upgrades,  and  even 
integration  of  future  systems. 

Defence  Research  Establishment  Ottawa  (DREO)  is  engaged  in  a  project  with  the 
Communications  Research  Centre  (CRC)  to  develop  a  proof-of-concept  testbed  for  a 
reconfigurable  omniband  (ROBR)  satcom  groimd  terminal  using  the  software  radio  concept 
described  above.  The  hardware  may  consist  of  general  processors,  digital  signal  processor 
(DSP)  boards,  application  specific  integrated  circuits  (ASICs),  or  field  programmable  gate 
arrays  (FPGAs)  to  perform  the  processing  functions  of  the  particular  waveform  standard  of 
interest. 

One  common  function  that  appears  in  the  transmit/receive  chain  of  a  groimd  terminal  is 
filtering.  Filter  specifications  may  differ  from  location  to  location  in  the  transmit/receive 
chain.  It  would  be  useful  to  have  one  software  module  capable  of  generating  filter 
coefficients  for  many  types  of  filters.  While  many  commercial  software  packages  are 
currently  available  for  filter  design,  they  produce  a  text  file  with  coefficients  that  have  to  be 
manually  integrated  into  the  processing  elements  of  the  terminal.  Any  change  in  filter 
specifications  would  require  a  new  text  file  to  be  generated  and  integrated.  The  development 
of  a  reconfigurable  digital  filter  module  for  the  ROBR  project  allows  the  coefficients  to  be 
generated  or  updated  while  the  terminal  continues  to  operate. 


1.1  Scope 

The  work  documented  in  this  report  provides  a  suite  of  several  kinds  of  reconfigurable  low- 
pass  filters  for  use  in  a  DSP  and  more  specifically,  in  the  ROBR  terminal  testbed.  The 
reconfigurable  digital  filters  have  been  programmed  in  C.  Infinite  impulse  response  (HR)  and 
finite  impulse  response  (FIR)  filters  have  been  implemented  for  more  flexibility.  The 
techniques  to  design  HR  filters  are  based  on  conversion  of  continuous-time  IIR  systems  into 
discrete-time  HR  systems.  This  project  explores  and  implements  two  different  conversion 
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methods:  the  Bilinear  Transformation  and  the  Impulse  Invariance  method.  Five  major  types 
of  IIR  filters  have  been  implemented:  the  Butterworth  filter,  the  Chebyshev  filter,  the  inverse 
Chebyshev  filter,  the  Elliptical  filter  and  the  Bessel  filter.  In  contrast  FIR  filters  are  almost 
entirely  restricted  to  discrete-time  implementations.  FIR  filter  design  is  an  approximation  of 
an  ideal  frequency  response  using  specific  approximation  methods.  Three  approaches  have 
been  explored  and  implemented  to  design  optimal  equiripple  filters.  They  are  the  “Frequency 
Response  Sampling  Design”  method;  the  “Design  by  Windowing”  method;  and  the  “Parks- 
McClellan”  algorithm.  The  first  two  methods  consist  of  sampling  the  frequency  response  and 
performing  an  inverse  discrete  Fourier  Transform  (IDFT)  to  compute  the  filter’s  coefficients. 
The  Parks-McClellan  method  uses  techniques  from  the  approximation  theory.  The  following 
report  provides  a  description  of  IIR  and  FIR  systems  and  filters.  The  report  describes  the 
various  filter  design  methods  mentioned  above  and  their  implementation  for  a  general 
purpose  processor.  The  report  also  describes  the  adaptation  of  two  of  the  implemented 
methods  for  use  on  a  DSP.  A  comparison  of  the  generated  filter  coefficients  for  each  filter 
with  MATLAB  implementations  is  presented.  The  report  also  discusses  implementation 
issues  related  to  digital  filter  design.  The  reconfigurable  digital  filter  modules  are  available  in 
the  package  digital Jilters.  Further  implementation  details  can  be  found  in  the  user’s  guide 
in  the  digitaljilters  package. 
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2  Digital  Filters 


Filters  are  an  important  class  of  systems  in  signal  analysis,  signal  processing  and 
communication.  A  filter  can  be  described  as  a  discrete-time  system  or  an  analog  system  that 
passes  certain  frequency  components  while  rejecting  others.  In  a  more  general  context,  any 
system  that  modifies  certain  frequencies  relative  to  others  is  also  called  a  filter.  The 
following  sections  present  the  fundamental  concepts  involved  in  discrete-time  or  digital 
filtering. 

2.1  Discrete-time  signal 

A  discrete-time  signal  is  an  indexed  sequence  of  real  or  complex  mraibers  denoted  by  x[n] 
[1].  It  is  a  function  of  an  integer-valued  variable,  n,  that  represents  an  instant  in  time.  In 
practice,  discrete-time  signals  are  derived  by  sampling  a  continuous-time  signal  Xc(t)  to 
produce  a  sequence  of  samples.  Alternatively,  a  sequence,  x[n],  can  be  represented  as  a  sum 
of  scaled,  delayed  impulses  [1]  as  follows 


x[w]= 

irss-flo 

where 


(1) 

(2) 


1 

0 


w  =  0 
n^O 


(3) 


Figure  1  shows  an  example  of  a  discrete-time  signal,  x[n],  representing  an  arbitrary 
continuous-time  signal,  Xc(t). 


Figure  1  Graphical  representation  of  a  discrete-time  signal 
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2.2  Linear  Time-Invariant  (LTI)  systems 

This  section  introduces  LTI  systems,  which  are  very  important  for  rmderstanding  discrete 
filtering.  LTI  systems  may  be  described  in  terms  of  the  effect  they  have  on  discrete-time 
signals.  Figure  2  shows  a  block  diagram  of  an  LTI  system.  The  input  x[«]and  the  output 
^[njof  the  LTI  system  are  discrete-time  signals.  An  LTI  system  may  be  viewed  as  a  black 
box  where  its  output  is  related  to  its  input  by  the  impulse  response  (or  transfer  fimction)  hk[n] 
of  the  system. 


Figure  2  Representation  of  a  Linear  Time-Invariant  System 

The  impulse  response,  hk[n],  is  the  response  of  the  system  to  an  impulse  8[n-k].  An  LTI 
system  can  be  completely  characterized  by  its  impulse  response.  We  can  obtain  the  output  of 
the  system  from  any  input  by  computing 

yb]= 

which  is  commonly  called  a  convolution  sum  [1]  and  is  denoted  by 

y[n]  =  x[n]*h^[n]  (5) 

2.3  The  Laplace  transform  and  the  Z-transform 

The  Laplace  transform  is  one  of  the  most  important  tools  used  in  si^ial  analysis.  It  is  used  to 
represent  the  frequency  spectrum  of  a  given  signal.  In  fact,  the  Fourier  transform  is  a  special 

case  of  the  Laplace  transform.  The  Laplace  transform  of  a  function,  x(/),  defined  for 

t  €  [-  oo,oo],  is  given  by 

00  * 

Z(5)=  jx(t)e-^'dt  (6) 

-oo 

where  j  is  a  complex  frequency  variable,  s  =  <t  +  y'Q  [2].  If  x(t)  describes  the  behaviour  of 
a  system  in  the  time  domain,  X(s)  represents  the  behavior  of  the  same  system  in  the 
complex  frequency  domain.  It  is  noted  that  for  the  purposes  of  this  report,  the  symbol  Q  is 
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used  to  represent  the  continuous  frequency  variable  whereas  the  symbol  o)  is  used  to  denote 
the  discrete  frequency  variable  as  will  be  described  further. 

The  equivalent  transformation  for  a  discrete-time  system  is  the  Z-transform.  The  Z-transform 
changes  the  representation  of  a  discrete  signal  from  the  time  domain  to  the  discrete  frequency 
domain.  The  Z-transform  of  a  discrete  signal  x[«]  is  given  by  [2] 


^(z)= 

In  order  to  see  the  relationship  between  the  Laplace  and  Z-  transforms,  consider  a  function, 
which  is  obtained  by  sampling  a  continuous  function,  represented 

mathematically  as 


-  nT)  =  'ZxAnTm  -xT)  (8) 

«»-oo 


where  T  is  the  sampling  period.  The  Laplace  transform  of  Xg(f)  can  be  written  as 


-y,W= 


(») 


»=-«> 


(10) 


Using  the  relationship  in  Equation  (2),  the  Z-transform,  X(z),  can  be  compared  with  the 
Laplace  transform,  (js) ,  where  it  can  be  seen  that  they  are  related  by  a  variable  change 

z  =  (11) 


so  that 

=  jf.w  ( 12) 


It  is  noted  that  the  substitution  z  =  transforms  the  5  =  JO,  axis  of  the  complex  frequency 
plane  onto  a  unit  circle  z  =  e-'™  [2]. 
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2.4  Filter  design  and  module  specifications 

Filters  are  an  important  class  of  LTI  systems.  They  are  used  extensively  in  communications 
(e.g.  low-pass  filters).  It  is  convenient  to  characterize  a  filter  by  its  frequency  response 
expressed  by  the  magnitude  of  its  transfer  function,  .  An  example  of  the  magmtude 

response  of  the  transfer  function  for  a  low-pass  filter  is  shown  in  Figure  3.  Parameters  that 
describe  the  filter  characteristics  are  listed  in  Table  1. 


Figure  3  Magnitude  response  for  a  low-pass  filter 


Filter  order 

n 

Passband  frequency 

Stopband  frequency 

a 

Passband  ripple 

6p 

Stopband  ripple 

8s 

Cutoff  frequency 

Oc 

Sampling  frequency 

Fs 

Passband  attenuation 

Ap 

Stopband  attenuation 

As 

Transition  band 

^s] 

Table  1  List  of  specifications  for  a  low-pass  filter 


Mathematically,  the  low-pass  filter  in  Figure  3  can  be  described  [1]  by 
\-d^<\H[e^^\<\  +  S^,  0<|Q|<i:2p 


(13) 
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(14) 


The  design  process  of  a  filter  begins  with  the  filter  specifications,  which  include  the 
constraint  on  the  magnitude  of  the  frequency  response,  the  type  of  filter  and  the  filter  order. 
Once  the  specifications  have  been  defined,  the  next  step  is  to  find  a  set  of  filter  coefficients. 
The  coefficients  are  simply  the  values  taken  from  the  transfer  function  AJn]  for  specific 

indices  that  produce  the  acceptable  filter  response.  After  the  coefficients  have  been 
generated,  the  next  step  is  to  use  them  to  process  a  signal. 

There  exist  two  classes  of  digital  filters:  infinite  impulse  response  (HR)  filters  and  finite 
impulse  response  (FIR)  filters.  They  are  both  described  in  more  detail  in  the  following 
sections. 
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3  Infinite  Impulse  Response  Filters 


3.1  Infinite  impulse  response  (HR)  LTI  systems 


HR  systems  are  a  subclass  of  linear  time  invariant  systems  and  satisfy  a  linear  constant- 
coefficient  difference  equation  [1]  of  the  form 

Y,a^yin-k)  =  J^b^x(,n-k)  Oq^I  (15) 

Applying  the  Z-transform  to  both  sides  of  Equation  (15)  and  using  the  linearity  and  the  time- 
shifting  properties  of  Z-transforms  [  1  ]  gives 


it«0  k^O 


The  transfer  function  in  the  Z-domain,  H(z),  can  now  be  defined  as 


Yiz) 

Z(z)  ,  ^ 

^  ^  1+2-"* 


.-k 


bf)  +6)7'*  +b2Z~^  +  ...  +  bi^z~^ 
l  +  a,z'‘  +a2z"^  +...  +  affZ~'^ 


(16) 


(17) 


The  coefficients,  6*  and  ,  are  the  filter  coefficients. 


3.1 .1  Flow  diagrams  of  recursive  structures 

For  implementation  on  a  general  or  digital  processor,  HR  systems  of  the  form  presented  in 
Equation  (15)  must  be  converted  into  a  structure  from  which  an  algorithm  can  be  derived. 
The  difference  equation  given  by  Equation  (15)  can  be  represented  graphically  as  a  recursive 
structure  [1].  The  HR  structure  shown  in  Figure  4  is  referred  to  as  Direct  Form  I. 


Figure  4  Flow  diagram  implementing  the  Direct  Form  I  realization  of  an  HR  filter 
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Table  2  shows  the  computational  requirements  of  Direct  Form  I  [3]. 


Number  of  multiplications 

N  +  M  +1  per  output  sample 

Number  of  additions 

N  +  M  per  output  sample 

Number  of  delays 

N  +  M 

Table  2  Computational  requirements  for  Direct  Form  I  structure 


3.2  HR  filter  design 

There  are  two  general  approaches  used  to  design  HR  filters.  The  first  approach  is  to  design  an 
analog  HR  filter  and  then  map  it  into  an  equivalent  digital  filter.  This  method  is 
computationally  efficient  and  gives  a  lot  of  control  on  the  design  because  the  art  of  analog 
filter  design  is  highly  advanced  [1].  The  second  approach  is  to  use  algorithmic  and  iterative 
design  procedures,  which  gener^ly  requires  solving  a  set  of  linear  or  non-linear  equations. 
The  first  approach  has  been  used  in  this  project. 

3.2.1  Analog  filter  design 

3.2.1. 1  Classical  analog  filter  approximations 

This  section  presents  the  characteristics  of  the  five  analog  filter  prototypes  that  have  been 
used  to  produce  discrete  HR  filters.  The  important  characteristics  to  be  considered  are  the 
transfer  function  of  the  filter,  the  magnitude  and  the  phase  of  the  fi'equency  response  of  the 
filter,  the  poles  and  zeroes  of  the  filters.  The  transfer  function  of  the  analog  prototypes  will 
be  expressed  in  the  continuous  frequency  domain  in  terms  of  its  Laplace  Transform. 

The  transfer  function  of  any  analog  prototype  filter  may  be  expressed  as 


-  (18) 


where  zi  are  the  zeroes  and  pk  are  the  poles  of  the  transfer  function.  The  magnitude  of  the 
fi:equency  response  may  be  expressed  as 


(19) 
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The  poles  of  the  transfer  function  [1]  are  complex  numbers  and  are  usually  of  the  form 


(20) 


3.2.1. 2  Butterworth  filter  properties 

A  Butterworth  filter  yields  a  flat  frequency  response  in  the  passband  and  in  the  stopband  as 
shown  in  Figure  5.  The  Butterworth  approximation,  generally  used  to  design  low-pass  filters, 
yields  an  allpole  filter  and  can  be  described  by  the  following  equations  [4]  [5]. 


normalized  transfer  fijnction 


unnormalized  transfer 
function 


magnitude 


H(s)  = 


n(5-p*) 

Jfc=l 


Q" 


ife-1 


2r\2n 


(21) 


(22) 


(23) 


/,  s  =  ^|: 


where  Qc  =  cutoff  frequency,  5  =  v  10  - 1 ,  and  Ap  =  passband  attenuation  in  dB. 

The  poles  of  H(s) ,  ,  are  located  at  2n  equally  spaced  points  around  a  circle  of 

radius  [4]  [5]  and  are  given  by 


A  =  (-l)‘^^"(y^^)  =  exp«{  j 


.(n  +  l  +  2k)7i\ 


In  J 


k  =  0,1,...,«-1 


(24) 


cr^  =  cos 


{n  +  \  +  2k)7t\ 

<  2n  J 


(25) 


/-V  •  I  {n  +  \  +  2k)7t 

Q,=Q.sin - - 


(26) 
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Figure  5  Frequency  response  for  a  Butterworth  filter 


The  number  of  poles  equals  the  order  of  the  analog  filter,  n.  As  shown  in  Figure  5,  the 
frequency  response  in  the  transition  band  becomes  steeper  as  the  filter  order  increases. 

The  minimum  order  that  will  ensure  an  attenuation  of  As  or  more  at  frequencies  Qs  and  above 
can  be  obtained  [5]  by  using 


log(lO'^^^'°  -l) 


f 

21og 

V 


Q 

Q 


.1 


cj 


where  Qs  =  stopband  frequency 


(27) 


3.2.1. 3  Chebyshev  filter  properties 

Chebyshev  filters  are  designed  to  have  an  amplitude  response  with  relatively  sharp  transition 
from  the  passband  to  the  stopband.  This  sharpness  is  achieved  at  the  expense  of  ripples  that 
are  introduced  into  the  response.  As  with  the  Butterworth  approximation,  the  Chebyshev 
approximation  yields  an  allpole  filter.  Figure  6  shows  examples  of  Chebyshev  filters  of 
various  order,  n.  As  the  order  increases,  the  ripple  in  the  passband  increases.  However  the 
tradeoff  is  a  sharper  transition  from  the  passband  to  the  stopband. 
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Figure  6  Frequency  response  for  a  Chebyshev  filter 


The  transfer  function  and  magnitude  response  of  a  Chebyshev  filter  [4]  [5]  is  given  by 


Transfer  function 
Magnitude  response 


His)  =  H, 


n 


*=0 


-Pk 

i^-pk) 


1 

i+sXH^) 


(28) 

(29) 


where  the  static  gain.  Ho,  is  given  by 


H. 


A=1 


t=i 


nodd 

neven 


(30) 


The  parameter,  s,  is  dependent  on  the  passband  ripple,  =  201og(<5p) ,  as  follows 


and  Tn(Q)  is  the  Chebyshev  polynomial  given  by 


cos(«cos‘‘(Q)) 

cosh(«cosh“'(fi)) 


0<Q<1 

Q>1 


(32) 


The  2n  poles  of  a  Chebyshev  filter  [4]  [5]  are  given  by  Equation  (29)  where 

k  =  0,1,."»«-1 


(Ji;  =  Qj,  cosi 


=Qj,sin| 


2n 


and 


r  = 


n 

2. 

2 

7C 

(i/r)+r 

y 

2 

1 - rV'” 

i+Vi+f 

k  =  1 


(33) 


(34) 


(35) 


3.2.1 .4  Inverse  Chebyshev  properties 

The  inverse  Chebyshev  approximation  yields  a  filter  which  has  zeroes  and  poles.  Depending 
on  whether  the  order  is  even  or  odd,  the  filter  will  have  as  many  zeroes  as  poles,  or  n  -  1 
zeroes  and  n  poles.  Figure  7  shows  the  frequency  response  of  an  inverse  Chebyshev  filter  for 
different  orders.  Just  as  the  Butterworth  and  Chebyshev  filters  showed,  the  frequency 
response  becomes  steeper  in  the  transition  band  as  the  order  of  the  filter  increases.  However, 
the  inverse  Chebyshev  filter  frequency  response  exhibits  a  flat  response  in  the  passband  and 
ripples  in  the  stopband.  The  ripples  in  the  stopband  increase  as  the  filter  order  increases. 
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Figure  7  Frequency  response  for  Inverse  Chebyshev  filters 


(36) 


(37) 


where  Tn(f^)  is  again  the  Chebyshev  polynomial  given  by  Equation  (32),  and  where  8 
depends  on  the  stopband  ripple  in  dB,  =  201og(^j)  as  follows. 


The  transfer  function  of  an  Inverse  Chebyshev  filter  [4]  is  given  by 


ff(s)=fr2s.^zi_ 

(s-a^) 


at  =■ 


Sk 


and  its  magnitude  response  is 


1^(7^^)!'  = 


i+^^[7;^(Q)r 


(38) 
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The  poles  of  an  Inverse  Chebyshev  filter  [4]  are  given  by  Equation  (36)  where 


o‘*  = 


COS] 


(2/  +  1)  £ 

2«  "^2 


sinh// 


-1 


/  =  0,1,.">”~1 


— 


Q,  sin 


(2/  +  1)  n 
In 


cosh  // j  i  =  0,1,—. «  - 1 


where 


M  = 


sinh''(l/g~') 


n 


(39) 


(40) 


(41) 


The  zeroes  of  an  Inverse  Chebyshev  filter  [4]  are  given  by 
For  n  odd 


For  n  even 


j 

(2A:  +  l);r^ 
In  J 


k  = 
k  ^  n 


/ 

cos 

V 


j 

Clk  +  \)7t'^ 
2n  J 


k  =  0,l,...,n-l 


(42) 


(43) 


3.2.1. 5  Elliptical  filter  properties 

By  allowing  ripples  in  the  passband,  Chebyshev  filters  obtain  better  frequency  selectivity 
than  Butterworth  filters  because  of  the  sharper  transition  band.  Elliptical  filters  improve  upon 
the  performance  by  permitting  ripples  in  both  the  passband  and  stopband.  Figure  8  shows  the 
frequency  response  of  an  Elliptical  filter  for  different  values  of  filter  order. 
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Figure  8  Frequency  response  for  an  Elliptical  Alter 
The  transfer  function  of  an  Elliptical  filter  [4]  [5]  is  given  by 


where 


His)  = 


-^0  TT  +  ^0< 

^1/'^  ^0; 


(44) 


fn-l 
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for  n  odd 
for  n  even 


(45) 


=  (46) 

0 [1  for  n  even  '  ' 

and  Off,  is  a  constant  defined  in  Appendix  B. 

The  analog  static  gain,  Ho,  can  be  calculated  [4]  using 
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‘^.n- 

for  n  odd 

TT 

for  n  even 

irflo, 

(47) 


and  the  magnitude  response  [4]  [5]  using 

1 

l  +  (48) 

R„{Q)  is  the  Chebyshev  Rational  Function  with  respect  to  the  centre  frequency 
=  Qq  =  1  and  where  s  is  the  ripple  factor  which  depends  on  the  passband  ripple,  5p, 
or  on  the  stopband  ripple  in  decibels,  rj  [3]  as  follows 


e  = 


I 


(49) 


The  Rational  Normalized  Function  R„(Q)  is  given  by  [2] 


R„(Q)  = 


(«-lV2  Q  2  _q2 


Jbr n  odd 


for  n  even 


(50) 


The  Elliptical  filter  poles  can  be  computed  following  a  calculation  algorithm.  Steps  to 
calculate  the  values  and  of  the  transfer  function  are  taken  from  [4]  and  are  listed 

in  Appendix  B.  Once  these  values  are  calculated  it  is  easy  to  obtain  the  poles,  the  zeroes  and 
the  gain  of  the  analog  filter  which  are  then  used  to  compute  the  coefficients  of  the  filter  [4]. 
The  important  parameters  that  have  to  be  considered  for  the  design  of  an  analog  Elliptical 
filter  are 


Qp  =  passband  frequency 
=  stopband  fi:equency 
Ap  =  maximum  passband  loss  (dB) 

=  maximum  stopband  loss  (dB) 
k  =  selectivity  factor  =0^/0^ 

Using  the  quadratic  formula,  the  i**'  pair  of  complex  pole  values  [4]  can  be  expressed  as 
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(51) 


The  zeroes  occur  at 


Pi  = 


-46o, 


2 


(52) 


3.2.1. 6  Bessel  filter  properties 

Bessel  filters  are  designed  to  have  maximally  flat  group-delay  characteristics.  As  a 
consequence,  there  is  no  ripple  in  the  impulse  and  step  response.  On  the  other  hand,  the 
rolloff  of  the  frequency  response  is  more  gradual,  making  the  transition  band  wider  than  for 
Butterworth  or  Chebyshev  type  filters.  The  Bessel  filter  is  also  an  allpole  filter.  Figure  9 
shows  the  frequency  response  of  Bessel  filters  for  various  values  of  n. 


Figure  9  Frequency  response  for  a  Bessel  filter  for  different  values  of  filter  order,  n 
The  transfer  function  of  a  Bessel  filter  [4]  is  expressed  as  follows 

Zi'W  =  -TT  (53) 


19 


where 


(54) 


^  2{n-k)\ 

*  “  2"-'‘k\in-k)\ 


(55) 


The  following  recursion  formula  [5]  is  used  to  determine  q„{,s)  from  and  q„.2(s) 


=(2«-l)^„_, 


(56) 


Using  this  recursion  formula  involves  the  computation  of  the  poles  of  the  analog  transfer 
function  with  a  numerical  analysis  formula.  These  poles  are  the  roots  of  Equation  (53),  Using 
these  formulas  involves  more  complexity  in  the  algorithm  and  there  is  no  guarantee  that  the 
computation  of  the  roots  will  converge.  A  very  efficient  way  to  implement  the  Bessel  filter  is 
by  using  pre-computed  values  for  its  poles.  The  Bessel  filter  was  implemented  using  a  look 
up  table  for  filter  orders  of  up  to  n  =  25.  Poles  of  the  analog  transfer  function  of  the  Bessel 
filter  have  been  computed  using  the  Signal  Processing  toolbox  from  MATLAB  and  are  used 
in  this  project.  The  analog  poles  are  converted  to  discrete  values  using  one  of  the  analog-to- 
digital  conversion  methods  shown  in  the  next  section  to  obtain  a  discrete  design  for  the 
implementation. 


3.2.2  Conversion  of  analog  HR  filters  for  digital  implementation 

This  section  describes  two  analog-to-digital  conversion  methods  that  have  been  used  for  HR 
filter  implementation  in  this  project.  The  idea  of  the  conversion  methods  is  to  transform  the 
analog  transfer  function  expressed  in  the  analog  frequency  domain  or  the  s-plane  into  a 
discrete  transfer  function  expressed  in  the  discrete  frequency  domain  or  the  z-plane. 

3.2.2.1  Bilinear  transformation 

The  bilinear  transform  [1]  is  a  mapping  from  the  s-plane  to  the  z-plane  defined  by 


H{z)  =  H, 


*2 

T 

where  Hc(s)  is  the  Laplace  transform  of  a  continuous  function. 


(57) 
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In  the  previous  sections,  the  equations  for  the  analog  poles  of  the  transfer  functions  of 
different  UR  filters  expressed  in  Ae  S-domain  were  presented.  To  compute  the  discrete  poles 
of  the  digital  HR  filter  from  the  analog  poles  of  the  analog  prototype  the  following 
relationship  is  used. 


21-z'’ 

5  = - 7 

n+z-’ 


(58) 


or  conversely, 

2-Ts 
^~2  +  Ts 


(59) 


where  T  is  the  sampling  period. 


3.2.2.1.1  Frequency  warping  function 

In  bilinear  transformation,  the  relationship  between  the  analog  or  continuous  frequency, 
and  the  digital  frequency,  co,  in  the  Z-domain  [1]  is  given  by 

n  =  |tan(|)  (60) 

Conversely,  to  compute  the  digital  frequency  corresponding  to  an  analog  frequency  in  the  S- 
domain,  the  following  equation  applies 


QT 


fi)  =  2arctaiJ 

\  2 


(61) 


Equation  (60)  is  plotted  and  shown  in  Figure  10.  It  is  noted  that  the  range  of  analog 
frequencies  -oo  ^  <<x3  maps  to  normalized  digital  frequencies  -n<  a)  The  frequency  o) 

=  7t  corresponds  to  half  the  normalized  sampling  frequency  F,  =  27t.  If  viewed  from  the  point 
of  view  of  Nyquist’s  Theorem,  since  the  bilinear  transformation  is  able  to  map  all  continuous 
frequencies,  Q,  to  digital  frequencies  below  FJ2,  there  are  no  aliased  components  resulting 
from  the  conversion  of  an  analog  filter  to  a  digital  one.  However,  the  effect  of  this 
transformation  is  the  non-linear  compression  of  the  discrete  frequency  axis  as  shown  in 
Figure  10. 
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Figure  10  Mapping  of  continuous  frequency  to  digital  frequency 


Z.2.2.2  Impulse  invariance  method 

The  “Impulse  Invariance”  method  consists  of  sampling  the  impulse  response  of  a  continuous¬ 
time  system,  hc(t)  to  get  a  discrete  impulse  response,  h[n],  as  follows 


h[n]  =  Th,(nT)  (62) 

From  Nyquist’s  sampling  theory,  it  can  be  shown  [1]  that  the  frequency  response  of  the 
discrete-time  filter  is  related  to  the  frequency  response  of  the  continuous-time  filter  by  [1] 

However,  in  the  case  of  an  ideal  bandlimited  continuovis-time  filter, 

|Q|^;r/r  (64) 


|u>|  <  ;r  ( 65) 

It  follows  from  Equation  (64)  that  the  discrete-time  and  continuous-time  frequencies  are 
related  linearly  by 

a)  =  QT  (66) 


/f,(yQ)  =  0, 

which  reduces  Equation  (63)  to 


J  T 

V  y 
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In  practice,  the  frequency  response  of  a  continuous-time  filter  is  not  bandlimited  as  described 
in  Equation  (65).  As  a  result,  aliasing  can  occur  between  successive  terms  of  Equation  (63). 
The  presence  of  aliasing  has  an  effect  on  the  design  of  a  discrete-time  filter  using  the  Impulse 
Invariance  method.  However,  if  the  fi’equency  response  of  the  continuous-time  filter  at 
higher  fi^quencies  is  sufficiently  low  and  the  sampling  fi’equency  is  sufficiently  high  enough, 
then  aliasing  is  minimized. 

It  is  shown  in  [1]  that  the  transformation  fi'om  continuous-time  to  discrete-time  can  be 
achieved  by  considering  the  continuous-time  fi:oquency  response  as  a  partial  fiiaction  so  that 

(67) 

Atmay  be  obtained  [5]  by  computing, 

A -[(*-?, (««) 


Taking  the  inverse  Laplace  transform,  the  impulse  response  of  the  continuous-time  filter 
becomes 


N 


[0 


t<0 


(69) 


The  impulse  response  of  the  discrete-time  filter,  obtained  by  sampling  /ip(0fi'om  Equation 
(62)  becomes 


A[n]  =  n^T)  =  rf;  A,  (,e“^ )"«[«]  ( 70) 

k=l 

where  u[n]  is  the  unit  step  fimction. 

By  taking  the  Z-transform  of  Equation  (70),  the  discrete  transfer  function  of  the  filter  is 
therefore  given  by 


H(z) 


N 

=lTr 

*=i  i 


(71) 
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3.3  HR  filtering  module  implementation 
3.3.1  System 

The  HR  filters  implemented  consist  of  two  important  components,  the  coefficients 
computation  module  and  the  signal  processing  module.  The  coefficients  computation  module 
is  responsible  for  designing  the  filter  and  computing  the  coefficients.  In  order  to  verify  the 
coefficients  computation  module,  a  signal  processing  module  is  implemented  that  uses  the 
coefficients  generated  by  the  computation  module  to  filter  a  test  signal. 

The  coefficients  computation  module  operates  in  three  basic  steps  as  shown  in  Figure  11. 
First,  it  computes  the  analog  poles,  zeroes  and  static  gain  of  the  analog  design.  Secondly,  it 
transforms  the  analog  filter  design  to  a  digital  design  and  obtains  the  digital  poles,  zeroes  and 
static  gain.  The  final  step  consists  of  factorizing  the  digital  poles  and  zeroes  to  obtain  the 
filter  coefficients.  Two  analog-to-digital  conversion  methods  have  been  implemented  and  are 
compared. 


Figure  11  Block  diagram  of  the  implementation  of  the  reconfigurable  HR 
filtering  modules 

Each  of  the  steps  outlined  above  makes  use  of  a  header  file  called  “Filter .h”.  The  header  file 
contains  routines  for  complex  number  mampulation  and  polynomial  mampulation  and  is 
listed  in  Appendix  C.  Complex  munber  manipulations  include  computing  the  norm, 
conversion  from  Euler  representation  to  trigonometric  representation  using  the  formula, 
=cos®  +  ysine;,  multiplication  and  division  of  two  complex  numbers,  computing  the 
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square  root  of  a  complex  number,  and  computing  a  number  raised  to  a  complex  power,  . 
Polynomial  manipulations  include  the  multiplication  of  two  polynomials  and  computing  an 
expansion  formula  to  obtain  the  roots  of  a  polynomial. 

The  following  two  sections  describe  the  coefficients  computation  and  signal  processing 
modules  in  more  detail.  The  software  developed  for  this  project  is  assembled  in  a 
digital  Jilters  package  for  future  use  in  the  ROBR. 

3.3.2  HR  coefficients  computation  modules 

The  principle  function  of  these  modules  is  to  compute  the  constants  a,,  and  for  an  HR 

filter  given  by  Equation  (17).  Five  HR  coefficients  computation  modules  have  been 
implemented  based  on  the  analog  designs  presented  in  Section  3.2.  Table  3  lists  the 
implemented  modules. 


Filter  type 

Module  name 

Butterworth  filter 

butter 

Chebyshev  filter 

chebl 

Inverse  Chebyshev  filter 

cheb2 

Elliptical  filter 

Ellip 

Bessel  filter 

Bessel 

Table  3  HR  Filter  types  implemented  by  the  corresponding  modules 

More  information  about  how  to  use  each  module  can  be  found  in  a  user’s  guide  included  with 
the  software  developed.  All  modules  require  two  input  parameters,  the  digital  cutoff 
frequency  in  Hertz  and  the  sampling  frequency  also  in  Hertz  .  To  enable  the  analog  design 
computation  block  to  calculate  the  analog  values  of  the  analog  design  the  input  cutoff 
frequency  needs  to  be  normalized  and  mapped  as  in  Equation  (60).  Figure  12  shows  the 
implementation  of  these  operations. 


//  normalization  of  the  digital  cutoff  frequency  over  0  to  2n 
float  Wd  =  fc  *  PI  /  (Fs  /  2.0); 

//  computation  of  the  analog  frequencies 

float  Wc  =  2.0/{1.0/Fs)  *  tan{Wd/2.0); 


Figure  12  Cutoff  frequency  normalization  and  analog  frequency  computation 

Wd  is  the  normalized  digital  cutoff  frequency  received  by  the  modules.  Wc  is  the  analog 
cutoff  frequency  required  to  compute  the  analog  design,  and  Fs  is  the  sampling  frequency. 
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3.3.2.1  Analog  design  computation 

The  first  step  of  the  reconfigurable  HR  module  is  to  compute  the  analog  poles,  zeroes  and 
static  gain  for  a  desired  analog  filter  design.  Figure  13  shows  the  implementation  of  the 
analog  design  computation  block  for  the  Chebyshev  module  (chebl).  This  block  implements 
the  equations  shown  in  Section  3.2. 1.3.  The  variables,  p,  z  and  /T  are  the  arrays  that  contain 
the  poles,  the  zeroes  and  the  analog  static  gain,  respectively.  The  poles  are  computed  using 
Equations  (33)  and  (34)  and  the  analog  static  gain  K  (introduced  as  //o  in  Section  3.2.1)  is 
obtained  by  computing  Equation  (30).  Recall  from  Section  3.2.1  that  the  Chebyshev  filter  is 
an  allpole  filter  and  thus,  does  not  have  zeroes.  This  is  why  the  elements  of  array  z  are  set  to 
zero.  The  same  procedure  is  carried  out  for  the  other  IIR  filters  with  the  appropriate 
equations  for  poles,  zeroes,  and  static  gain. 

As  explain  earlier  in  Section  3. 2. 1.6,  the  Bessel  filter  is  implemented  using  a  lookup  table. 
No  zeroes  are  computed  because  this  filter  is  an  allpole  filter.  Figure  14  shows  the 
implementation  of  the  Bessel  coefficients  as  a  look  up  table.  The  module  considers  filters  up 
to  order  25.  The  poles  stored  in  the  table  have  been  calculated  with  MATLAB’s  besselap 
function  from  the  Signal  Processing  toolbox. 


//computation  of  the  analog  poles 
for  (i=l;i<=N;i++) 

{ 

z[i-l] .Q  =  0; 
z[i-l] .1  =  0; 

p[i-l].Q  =  sin ( {2*i-l) *PI/ (2*N) )*(( (1.0/gamma) -gamma) /2.0) ; 
//formula  for  analog  poles 

p[i-l).I  =  cos ((2*i-l)*PI/(2*N) )*(( (1.0/gamma) +gamma)/2.0); 
p[i-l]  =  set_to_zero(p[i-l] ) ; 

K  =  multc(K.Q,K.I,-l*p[i-l] .Q,-l*p[i-l] .1) ; 

//computation  of  the  analog  static  gain 

} 

//if  analog  static  gain  the  order  is  even,  adjust  the  analog  gain 
if (N%2==0) 

K.Q  /=sqrt((1.0  +  epsilon*epsilon) )  ; 


Figure  13  Chebyshev  analog  design  computation  algorithm 
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Switch (N) 

{ 

case  0:  return; 
case  1:  p[0].Q  =  -1.0; 

p[0] .1  =  0,0; 

break; 

case  2:p[0] .Q=-0 . 8660254037844386467637229; 

p[0] .I=+0. 4999999999999999999999996; 
p[l] .Q=-0. 8660254037844386467637229; 
p[l] .I=-0. 4999999999999999999999996; 
break; 

case  25:  p[0].Q=0.0; 

p[0] .I=-0.9062073871811708652496104; 
p[l] .Q=-0. 9028833390228020537142561; 
p[l] .I=-93077131185102967450643820.0E-27; 
p[2] .Q=-0. 9028833390228020537142561; 

} 


Figure  14  Implementation  of  the  Bessel  analog  design  computation  block  as  a  lookup 
table 

3.3.2.2  Analog-to-digital  Conversion 

The  analog-to-digital  conversion  block  is  the  2"**  step  in  the  coefficients  computation  module. 
Figure  15  shows  the  implementation  of  the  conversion  block  code  based  on  the  Bilinear 
Transform  method.  This  portion  of  code  belongs  to  the  function  “vo/V?  bilinearO”,  which  can 
be  found  in  the  file  “bilinear.c”  in  the  digital Jilters  package.  The  bilinear  conversion 
function  takes  the  analog  poles,  zeroes  and  static  gain  computed  in  the  previous  step  and 
computes  the  corresponding  digital  poles,  zeroes  and  static  gain. 

3.3.2.3  Factorization 

Factorization  of  the  zeroes  and  poles  of  the  discrete  transfer  function  obtained  from  the 
analog  design  is  done  using  an  expansion  formula.  Using  the  discrete-time  transfer  function 
given  in  Equation  (17),  the  poles  and  zeroes  of  the  HR  filter  are  the  roots  of  the  denominator 
and  numerator  polynomials  respectively.  The  expansion  formula  may  be  expressed  as 
follows 

^n+I  ~  ^/i+l  ^ 

where  c„  are  the  computed  coefficients  and  are  the  roots  of  the  polynomial.  The  variable 
n  is  the  order  of  the  polynomial  and  is  decremented  by  one  at  each  iteration.  The  variable  m 
equals  the  order  of  the  polynomial  (i.e.  n  =  m  for  the  first  iteration)  and  stays  constant  for  a 
set  of  iterations.  The  factorization  will  yield  n+1  coefficients  where  n  is  the  filter  order. 
Figure  16  shows  the  implementation  of  the  expansion  formula.  This  portion  of  code  belongs 
to  the  fimction  “voiVf  coeffO”  included  in  the  header  file  “Filter.h”. 
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//  discretization  of  the  zeroes 


for(i  =  0;i  <  M;i++) 
zd [ i ] =z [ i ] ; 

for(i  =  M;i  <  N;i++) 

{ 

zd [ i ] . Q=- 1.0; 
zd [i] . 1=0 . 0; 

) 

for(i  =  0;i  <  M;  i++) 

{ 

zd[i]  =  divc(2.0  +  (zd[i].Q  *  Wc/Fs  ),zd[i].I  *  Wc/Fs  ,2.0  -  zd[i] .Q 
*  Wc/Fs  ,-1.0  *  zd[i].I  *  Wc/Fs); 
zd[i]  =  set_to_zero(zd[i] ) ; 

tempNKd  =  multc (tempNKd.Q, tempNKd. I, (2 . 0*Fs/ (Wc) -z [i] .Q) , 

(-1.0*z[i] .1) ); 
tempNKd  =  set_to_zero (tempNKd) ; 

} 

//  discretization  of  the  poles 

for(i  =  0;i<  N;  i++) 

{ 

pd[i]  =  divc(2.0  +  (p[i).Q  *  Wc/Fs  ),p(i].I  *  Wc/Fs  , 

2.0  -  p[i].Q  *  Wc/Fs  ,-1.0  *  p(i].I  *  Wc/Fs); 
pd[i]  =  set_to_zero(pd[i] )  ; 

tempDKd  =  multc (tempDKd.Q, tempDKd. I,  (2*Fs/ (Wc) -p[i] .Q) , 


tempKd  =  dive (tempNKd.Q, tempNKd. I, tempDKd. Q, tempDKd. I) ; 
Kd  =  (multc (tempKd.Q, tempKd. I, K.Q, K. I) ) .Q; 


Figure  15  Bilinear  transform  analog-to-digital  conversion  algorithm 


num  =  multc  (c  [n]  .Q,  c  [n]  .  I,  e  [m]  .Q,  e  [m]  .  I ) ; 
c[n+l].Q  =  c[n+l].Q  -  num.Q; 
c[n+l].I  =  c[n+l].I  -  num.I; 
if (n  >  0) 

return  coeff (&c [0] , e, n-l,m) ; 

else 

return  0 ; 

Figure  16  Expansion  formula  algorithm 
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The  function  “vo/V/  coeffQ"  is  called  by  the  coefficients  computation  modules  after  the 
analog-to-digital  conversion  step.  The  function  will  fill  the  empty  array  cD[n+l],  where  n  is 
the  order  of  the  HR  filter  with  the  desired  filter’s  coefficients.  More  information  about  this 
subroutine  is  included  in  the  user’s  guide. 


3.3.3  HR  signal  processing  module 

The  HR  signal  processing  module  has  been  implemented  based  on  the  Direct  Form  I 
implementation  of  the  difference  equation  presented  in  Section  3.1.  This  signal  processing 
module  can  be  found  in  the  package  digital  Jilter  under  the  name  of  IIRDFL 

The  signal  processing  module  computes  the  output  of  a  system  using  the  IBR.  filter 
coefficients  generated  by  the  coefficients  computation  module.  The  input  test  signal  used 
was  an  impulse  sequence. 

Figure  17  shows  the  implementation  of  the  HR  signal  processing  algorithm.  Details  on  how 
to  execute  the  HR  filter  modules  and  signal  processing  module  are  described  in  Section  4.3.4. 
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for  (i=0;  i<nuincoef ;  i++) 

coeffden[i]  =  dive (aO . Q, aO . I, coef fden [i] .Q, coef fden [i] .1); 

for  (i-0;  i<numsainples;  i++) 

{ 

sumnuin .  Q=0 . 0 ; 
sumnum. 1=0 . 0; 
sumden.Q=0 .0; 
sumden. 1=0 . 0; 

for  ( j=0;  j<nuincoef ;  j++) 

indx  =  i  -  j ; 
if(indx  <  0)  break; 
else 

{ 

sumnum. Q  +=  coeffnum[ j ] .Q  *  samples [indx] ; 
sumnum. I  +=  coef fnum[ j ] . I  *  samples [indx] ; 

} 

} 

for ( j=l; j<numcoef ; j++) 

{ 

indx  =  i  “  j ; 

if (indx  <  0)  break; 

else 

{ 

sumden. Q+=multc (coef fden [j ] ,Q, coef f den [j ] .1/ 

output [ j-1] .Q/ output [ j~l] .1) .Q; 
sumden. I+=multc( coef f den [j] .Q, coef fden [j ] .1, 

output [j-1] .Q, output [j-1] .1) . I; 

} 

} 


for ( j=numcoef-l; j>0; j — ) 

output [ j ] =output [ j - 1 ] ; 

output [0] . Q=sumnum. Q-sumden . Q; 
output [0] . I=sumnum. I-sumden , I ; 
} 


Figure  17  IIR  signal  processing  module  algorithm 
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4  Finite  Impulse  Response  Filters 

4.1  FIR  LTI  systems 

The  FIR  signal  processing  module  implements  an  FIR  structure  called  the  Direct  Form  to 
process  digital  signals.  For  FIR  systems,  the  transfer  function  H{z)  has  no  poles  except  at 
z  =  0 .  Thus,  H{z)  is  simply  a  polynomial  in  z"'  of  the  form 

M 

H{z)  =  2;  ^*2"*  =  /%  +  + ...  + 

The  output  of  such  a  filter  is 

y(ri)  =  Y,Kk)x(n-k) 

=  /^,  +  1\x{n  - 1)  +  h^x{n  -  2)...h^x{n  -  M) 

which  is  the  computational  sum  introduced  in  Section  2.2. 

4.1 .1  Flow  diagrams  for  non-recursive  structures 

The  flow  diagram  shown  in  Figure  18  illustrates  the  convolution  sum  that  relates  an  FIR 
filter’s  output  to  its  input.  This  structure  is  called  the  Direct  Form.  Unit  delays  are  denoted  by 
z“*  as  shown  in  the  figure  below 


Figure  18  Flow  diagram  of  the  Direct  Form  realization  of  an  FIR  filter 


Table  4  shows  the  computational  requirements  of  Direct  Form  implementation  of  an  FIR 
filter  [3]. 


(73) 


(74) 
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Number  of  multiplications 

N  +  1  per  output  sample 

Number  of  additions 

N  per  output  sample 

Number  of  delays 

N 

Table  4  Computational  requirements  of  the  Direct  Form  I  structure  of  FIR  filters 


4.2  FIR  filter  design 

Recall  for  HR  filters  that  the  design  techniques  were  based  on  transformation  of  continuous¬ 
time  HR  systems  into  discrete-time  systems.  In  contrast,  FIR  filter  design  is  almost  entirely 
restricted  to  discrete-time  implementation.  Consequently,  the  design  techniques  for  FIR 
filters  are  based  on  directly  approximating  the  desired  frequency  response  of  the  discrete¬ 
time  system.  Furthermore,  most  of  these  approximation  techniques  avoid  the  problem  of 
factorization  that  complicates  the  design  of  HR  filters.  In  the  context  of  this  project,  three 
techniques  have  been  implemented  to  compute  the  FIR  filter  coefficients  of  interest:  the 
“Frequency  Sampling  Design”  technique;  the  “Design  by  Windowing”  method;  and  the 
“Parla-McClellan”  method.  The  following  sections  introduce  these  three  methods. 

In  addition,  the  design  of  a  Gaussian  filter  and  a  digital  integrator  is  included  in  this  section. 

A  Gaussian  filter  and  a  digital  integrator  are  used  in  the  premodulation  stage  of  a  Gaussian 
minimum  shift  keying  (GMSK)  modulator  that  is  to  be  implemented  in  the  ROBR. 

4.2.1  “Frequency  Sampling  Design”  filter  module 

The  “Frequency  Sampling  Design”  method  is  a  straightforward  design  procedure.  The 
frequency  response  of  an  ideal  filter  is  sampled  and  each  sample  of  the  frequency  response  is 
a  coefficient.  Figure  19  shows  the  frequency  response  of  an  ideal  low-pass  filter.  To  get  the 
values  of  the  coefficients  in  the  time  domain,  an  IDFT  is  performed  on  the  samples  collected 
[1][3][5]. 


Figure  19  The  desired  frequency  response  Hj{(o)  of  an  ideal  low-pass  filter 
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The  desired  frequency  response  is  uniformly  sampled  at  N  equally  spaced  points  between  0 
and  271  to  yield 


Hik)  =  H, 


f 


k  =  0X...,N-l 


(75) 


These  samples  constitute  an  N-point  DFT,  whose  inverse  is  the  impulse  response  of  an  FIR 
filter  of  order  N-1: 


The  inverse  DFT  can  be  modified  to  take  advantage  of  symmetry  conditions.  Table  5  shows 
adapted  IDFT  formulas  to  the  four  types  of  FIR  filters.  [5] 


Type 

hinl  for  n  =  0, 1, 2, ...,  N-1 

1 

h[n]  symmetric 
Nodd 

2 

h[n]  symmetric 
N  even 

3 

h[n]  asymmetric 
Nodd 

4 

h[n]  asymmetric 
N  even 

Table  5  Inverse  Discrete  Fourier  Transform  formulas  for  FIR  Design 

In  this  project,  the  first  type  and  the  second  type  of  FIR  filter  have  been  implemented  for  a 
low-pass  filter  with  synunetric  impulse  response. 


4.2.1. 1  Gibbs  phenomenon 

One  problem  related  to  the  approximation  methods  used  to  produce  FIR  filter  coefficients  is 
Gibbs  phenomenon.  Figure  20  shows  the  amplitude  of  the  frequency  response  and  the  power 


33 


spectrum  of  a  low-pass  FIR  filter  affected  by  Gibbs  phenomenon.  The  coefficients  of  the 
filter  were  computed  with  the  “Frequency  Sampling  Design”  method.  The  frequency 
response  has  an  oscillating  behaviour  that  is  more  pronounced  near  the  edge  of  the  passband. 
This  behavior  is  known  as  Gibbs  phenomenon  and  is  the  result  of  approximating  a 
discontinuity  in  the  desired  frequency  response.  In  [3]  it  is  noted  that  if  a  fimction  with  a 
discontinuity  is  approximated  by  a  Fourier  series,  there  is  an  overshoot  in  the  region  near  the 
discontinuity.  As  the  number  of  Fourier  series  terms  increases,  the  squared  error  decreases 
and  approaches  zero  as  the  number  of  terms  approaches  infinity.  However,  the  maximxun 
value  of  the  overshoot,  and  therefore  the  maximum  value  of  the  error,  do  not  go  to  zero  but 
approaches  a  constant  value  of  1 1  %  of  the  size  of  the  discontinuity.  In  the  “Frequency 
Sampling  Design”  method,  the  overshoot  may  approach  approximately  18%  of  the 
discontinuity  [3]. 


Figure  20  (a)  Frequency  response  of  an  FIR  filter  affected  by  Gibbs  phenomenon 

(b)  Power  Spectrum  of  an  FIR  filter  affected  by  Gibbs  phenomenon 


4.2.3  “Design  by  Windowing”  filter  module 

The  “Design  by  Windowing”  method  begins  with  a  desired  frequency  response  that  can  be 
represented  as 


««s-00 


(77) 


where  hj{n\  is  the  corresponding  impulse  response  sequence.  Let  be  an  ideal  low- 

pass  filter  with  fi’equency  response 
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(78) 


Q<n<M 

otherwise 


where  M  is  the  filter  order.  As  in  the  “Frequency  Sampling  Design”  method,  the  impulse 
response  of  the  ideal  frequency  response  can  be  obtained  by  performing  an  inverse  Fourier 
transform. 


hAn\  =  ^]HAe‘"")dm 


(79) 


However,  to  improve  the  impulse  response  of  the  filter  and  to  reduce  Gibbs  phenomenon  the 
ideal  impulse  response  is  truncated  using  a  window.  The  simplest  way  to  obtain  a  causal  FIR 
filter  from  [n]  is  to  define  a  new  system  with  impulse  response  h{n\  given  by 


O^n^M 

otherwise 


(80) 


where  M  is  the  order  of  the  transfer  function  polynomial.  Thus,  (M  + 1)  is  the  length  of  the 
impulse  response.  Alternatively,  we  can  represent  h[ri\  as  the  product  of  the  desired  impulse 
response  and  a  finite-duration  “window”,  >v[n] , 


h{n\  =  hj  [«]w[n] 


(81) 


This  multiplication  truncates  the  ideal  infinite  impulse  response,  hd[n],  to  obtain  a  finite 
impulse  response,  h[n],  with  less  imperfection,  thus  reducing  the  effects  of  Gibbs 
phenomenon.  In  the  frequency  domain.  Equation  (81)  can  be  expressed  as 


Hie^^^)  =  —  \H^{e^‘^)W(eJ‘^-^)dd  ( 82) 

That  is,  /f(e^‘“)is  the  periodic  convolution  of  the  desired  ideal  frequency  response  with  the 
Fourier  transform  of  the  window.  Thus,  the  frequency  response  will  be  a  “smeared” 

version  of  the  desired  response  Hj{e^‘“) . 

Some  commonly  used  windows  [1]  are  shown  in  Figure  21  and  their  equations  listed  in  Table 
6.  The  Hamming  window  used  for  this  project  is  shown  in  red  in  Figure  21 . 
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Commonly  used  windows 


Figure  21  Commonly  used  windows 


Window 

Window  equation 

Rectangular 

win)  =  ^ 

1,  0<n<M 

0  otherwise 

Barlett  (triangular) 

w{n)  = " 

'intM,  Q<n<MI2 

2-2nlM  M/2<n<M 

0  otherwise 

Hanning 

0.5-0.5cos(2m/ M)  0<n<M 

[O  otherwise 

Hamming 

w{n)  =  > 

[0.54- 0.46  cos(2;zn/M)  0<n^M 

[O  otherwise 

Blackman 

win)  =  > 

[0.42-0.5  cos(2;n7  / M)  +  0.08 cos(4;277  /  M)  0  ^  <  M 
[O  otherwise 

Table  6  Commonly  used  window  functions 

Figure  22  shows  the  improvement  obtained  with  the  “Design  by  Windowing”  method  as 
compared  with  the  “Frequency  Sampling  Design”  method.  The  overshoot  near  the 
discontinuity,  in  the  passband  and  in  the  stopband,  has  been  considerably  reduced  using  the 
“Design  by  Windowing”  method. 
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(a) 


Figure  22  Comparison  between  the  “Frequency  Sampling  Design”  and  “Design  by 
Windowing”  methods  for  FIR  filter  design,  (a)  Frequency  response,  (b)  Power 
spectrum 


The  power  spectrum  on  the  right  side  shows  that  the  attenuation  in  the  stopband  has  been 
improved  by  more  then  20  dB. 


4.2.4  Parks-McClellan  filter  module 

While  the  design  of  FIR  filters  with  the  “Frequency  Sampling  Design”  method  with  or 
without  windowing  is  straightforward,  there  are  a  number  of  limitations.  The  “Parks- 
McClellan”  algorithm  yields  optimal  filters  and  offers  more  control  on  certain  regions  of  the 
frequency  response  of  the  filter.  The  Parks-McClellan  algorithm  is  based  on  expressing  the 
filter  design  problem  as  a  problem  in  polynomial  approximation  [1]  which  is  described 
briefly  in  &e  next  section. 


4.2.4.1  Chebyshev  approximation 

Given  the  problem  of  designing  a  low-pass  filter  with  specifications  such  as  those  shown  in 
Figure  3,  the  Parks-McClellan  algorithm  considers  that  the  desired  j&equency  response  of  a 
filter  may  be  approximated  by  a  c’th-order  polynomial  in  cosioas  follows 

P(u>)  =  Ja^Ccosfi))*  (83) 

*=0 
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where  c  =  ,  N  is  the  number  of  samples  of  the  filter  impulse  response,  and  coefficients 

are  chosen  so  as  to  yield  a  P(<y)  which  is  optimal  in  a  sense  that  is  defined  below. 


Figure  23  The  desired  frequency  response  D((d)  of  a  low-pass  filter 


Consider  D(o))  to  be  the  desired  frequency  response,  also  shown  in  Figure  23,  so  that 


Dico)  = 


for  0<a}<o)p 
for  a)g<(o<n 


(84) 


and  let  W(co)he  the  weighting  function  for  the  approximation  error  over  each  of  the  intervals 
that  the  filter  is  defined. 


W(co)  = 


fl  for  0<CD<co^ 


6,1 


for  o)„<(o<7r 


(85) 


where  8i  and  82  are  the  amplitude  of  the  passband  and  stopband  ripple  respectively.  The  error 
made  by  the  approximation  of  D(co)  by  P(co)  can  be  computed  by  [1][6][7] 


E{a)  =  W{(D)[D{co)-P{(o)] 


(86) 


The  weighted  error  function,  E{o)),  the  weighting  function  W{a)),  the  desired  frequency 
response  D{a))  and  the  approximation  polynomial  P{co) ,  are  defined  for  the  same  discrete 
subset  of  frequencies  taken  from  the  interval  0  <  u;  <  ;r .  For  a  low-pass  filter,  these  four 
functions  will  be  defined  over  Q<(t)<C0p  and  (Og<(o<n  as  indicated  in  Equations  (84) 
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and  (85).  The  subinterval  is  the  transition  band  of  the  filter.  As  Eia>),  W{o))y 

i>(fi>)and  P{ct)),  are  not  defined  over  the  transition  band,  the  Parks-McClellan  method  [1] 

allows  P((o)  to  take  any  shape  in  this  band  to  achieve  its  optimum  approximation  and  meet 
the  filter  specifications.  The  criterion  used  in  this  design  procedure  is  the  Chebyshev 
Criterion  [1]  where  within  the  frequency  intervals  of  interest,  P(a))  is  chosen  to  minimize 
the  maximum  weighted  approximation  error.  This  can  be  expressed  by 


min  |maxl^(®)|j 

/i[ii]:0^^L  t  (oeF  J 


where  F  is  the  closed  subset  of  0  ^  <  ;r  over  which  the  filter  is  specified  (i.e.  0  ^  nf  ^  a) ^ 

and  <  o  ^  ;r ).  Thus,  the  method  is  to  seek  the  set  of  fi'equencies  that  produce  the  impulse 

response  values  that  minimize  the  error.  These  impulse  response  values  are  then  used  to 
compute  the  coefficients  of  a  desired  filter.  An  IDFT  will  be  performed  on  the  impulse 
response  values  to  give  the  coefficients  of  the  filter  in  the  time-domain.  A  weighting  factor 
may  be  associated  with  each  fi'equency  subinterval  (i.e.  0  to  and  oo^^oo-^ti  for  a 

low-pass  filter).  The  general  idea  for  using  a  weighting  factor  is  to  exaggerate  the  error  of  the 
approximation.  The  algorithm,  in  turn,  will  try  to  produce  a  better  approximation  for  bands 
with  the  higher  weighting  factor  by  minimizing  this  amplified  error.  This  results  in  more 
iterations  to  compute  the  approximation  but  a  more  accurate  response  is  achieved  in  the 
higher  weighted  band. 

Parks-McClellan  applies  the  Alternation  Theorem  from  approximation  theory  to  minimize 
the  error.  The  Alternation  Theorem  is  described  in  the  next  section. 


4.2.4.2  Alternation  theorem 

The  Alternation  Theorem  states  that  if  the  firequency  response  of  a  filter  is  represented  by  a 
linear  combination  of  c  cosine  functions,  as  expressed  in  Equation  (83),  the  best-weighted 
Chebyshev  approximation  to  a  desired  fi'equency  response,  Dijxi)  is  achieved  if  the  weighted 

error  function  E{(o)  exhibits  at  least  c  +  2  extremal  firequencies  inF^.  An  extremal 
firequency  is  the  fi'equency  at  which  a  ripple  (either  in  the  passband  or  the  stopband)  is  at  its 
Tnayimiim  or  minimum  value.  is  the  subset  of  frequency  intervals  over  which  the 

passband  and  stopband  are  defined.  Thus,  there  are  at  least  c  +  2  firequencies  in  F^  such  that 
6),  <  Oj  <  —  <  ®c+2  for  /  =  l,2,...,c  + 1  [1].  To  find  the 

extremal  fi'equencies  fi’om  the  discrete  subset  of  potential  extremal  fi*equencies  F^,  the 

Parks-McClellan  method  uses  the  Remez  Exchange  algorithm.  The  Remez  Exchange 
algorithm  is  a  set  of  conditional  statements  applied  to  the  potential  extremal  fi'equencies 
contained  in  F^ .  Based  on  the  Chebyshev  error  criterion,  the  Remez  Exchange  algorithm  will 

determine  if  a  fi’equency  is  an  extremal  frequency  or  not.  The  c+2  extremal  firequencies 
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found  by  the  algorithm  are  then  used  to  compute  the  filter  fi'equency  response  with  the 
approximation  fimctionP(fi)) .  The  Remez  Exchange  algorithm  is  described  further  in  [6]  [7]. 

From  the  Alternation  Theorem  we  can  write  [7] 

E{(a,)  =  W{co,)[D(,o)^)-P{,o),)]={-\t'S  ,  /  =  l,2,...,(c  +  2)  (88) 


where  S  is  the  optimum  error.  Parks  and  McClellan  [1][6][7]  found  that  for  the  given  set  of 
the  extremal  fi-equencies,  d  is  given  by  the  formula 


where 


and 


h 


c+2 


i.i  (•«.  -•«,) 

i^k 


=COSfi>, 


(89) 


(90) 


(91) 


Parks  and  McClellan  used  the  Lagrange  interpolation  formula  [1][6]  to  obtain 


where 


Xk/(^-^*)]c* 


X  =  cos  to 


C,=D{q})- 


W(o)) 


c+l 


^*=n 


1 


7-1  (J^i-Jc,) 


-  b^ipc^  Xi^y) 


(92) 


(93) 

(94) 

(95) 
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Figure  24  shows  an  example  of  a  low-pass  filter  described  by  the  Parks-McClellan  method. 
A  flow  diagram  of  the  algorithm  can  be  found  in  [1]. 


Figure  24  Typical  example  of  a  low-pass  filter  approximation  that  is  optimal  according 
to  the  alternation  theorem  for  c  =  7 


4.2.5  Gaussian  digital  filter  module 

A  digital  OMSK  modulator  is  to  be  implemented  in  the  ROBR  testbed.  One  implementation 
of  the  OMSK  modulator  requires  a  Gaussian  filter  and  a  digital  integrator  as  shown  in  Figure 
25.  As  a  result,  a  Gaussian  filter  is  also  implemented  using  an  FIR  filter  design  technique. 
The  coefficients  computation  module  uses  the  “Frequency  Sampling  Design”  method  to 
calculate  the  coefficients  of  a  Gaussian  filter.  However,  instead  of  sampling  the  ideal 
fi'equency  response  of  a  low-pass  filter,  a  Gaussian  distribution  is  sampled.  The  impulse 
response  of  a  Gaussian  filter  is  expressed  as  [12] 


5 


where  = 


n 

V21n2 


(96) 


or 
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cos(nw„) 


Figure  25  The  premodulator  stage  of  GMSK  digital  modulator  includes  a  digital 
integrator  followed  by  a  Gaussian  filter 

1.2 
1 

0.8 

H(C0)  0.6 

0.4 
0.2 
0 

normalized  frequency 

Figure  26  Frequency  response  for  a  Gaussian  filter 
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4.2.5.1  Digital  integrator  moduie 

As  mentioned  in  Section  4.2.5,  a  OMSK  modulator  is  to  be  implemented  in  the  ROBR  and 
requires  a  digital  integrator.  As  a  result,  the  implementation  of  a  digital  integrator  module  is 
included  here.  The  implementation  of  the  digital  integrator  is  very  simple.  The  output  for  a 
given  moment,  yk,  of  such  a  filter  is  the  summation  of  the  present  input  the 

previous  or  past  inputs,  which  can  be  expressed  mathematically  as 

k 

(9») 

i=0 

4.3  FIR  filter  module  implementation 

4.3.1  System 

As  with  the  HR  filter  modules,  the  FIR  filter  modules  implemented  consist  of  two 
components,  a  coefficients  computation  module  and  a  signal  processing  module.  The 
coefficients  computation  module  is  subsequently  broken  down  into  two  stages,  a  sampling 
step  and  the  computation  of  an  inverse  discrete  Fourier  Transform.  T^e  coefficients 
computation  module  takes  the  filter  specifications,  computes  the  desired  FIR  design  and 
generates  the  coefficients  of  the  filter.  The  signal  processing  module  takes  an  input  signal 
and  processes  it  using  the  coefficients  computed  by  the  coefficients  computation  module. 
Figure  27  shows  the  architecture  of  the  FIR  filter  modules. 


Input 

signal 


Output 

signal 


Figure  27  Block  diagram  of  FIR  filter  module  implementation 
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The  three  design  methods  described  earlier  have  been  implemented  to  compute  FIR  filter 
coefficients.  The  “Frequency  Sampling  Design”  method,  the  “Design  by  Windowing” 
method  and  the  “Parks-McClellan”  method  are  described  in  the  following  sections.  In 
addition,  the  implementation  of  filter  modules  for  ^e  Gaussian  filter  and  the  digital 
integrator  are  presented. 


4.3.2  FIR  coefficients  computation  modules 

For  FIR  filters,  the  coefficients  computation  modules  generate  the  coefficients  identified 
in  Equation  (73).  The  coefficients  are  subsequently  used  in  the  FIR  signal  processing 
modules  to  process  a  digital  signal.  Five  FIR  coefficient  computation  modules  have  been 
implemented.  The  implemented  modules  are  listed  in  Table  7. 


Filter  Design 

Module  Name 

Frequency  Sampling  Design  method 

Freqsampling 

Design  by  Windowing  method 

Hamming 

Park-McClellan  Method 

Remezex 

Gaussian  filter 

Gauss 

Digital  Integrator 

Dintegrator 

Table  7  FIR  filter  types  implement  by  the  corresponding  modules 


More  information  about  how  to  use  each  module  can  be  found  in  the  user’s  guide  included 
with  the  software  developed.  All  FIR  coefficient  computation  modules,  except  the  digital 
integrator  module,  require  two  input  parameters:  the  digital  cutoff  fi-equency  and  the 
sampling  frequency. 


4.3.2.1  Frequency  sampling  design  method  module  implementation 

The  design  method  starts  in  the  frequency  domain  and  is  separated  into  two  steps.  To  apply 
the  “Frequency  Sampling  Design”  method,  a  vector  is  created  which  contains  the  fiiequency 
response  samples  of  the  desired  low-pass  filter.  The  size  of  this  vector  will  be  the  same  as  the 
number  of  coefficients  to  be  computed.  This  operation  is  equivalent  to  sampling  an  ideal 
fi’equency  response.  Once  this  vector  is  built,  an  IDFT  is  then  performed  on  the  vector’s 
elements  to  compute  the  coefficients  of  the  filter  in  the  time  domain.  The  implementation  of 
these  two  steps  to  design  a  low-pass  filter  is  described  in  the  next  sections. 

4.3.2.1.1  Sampling  the  ideal  frequency  response 

The  first  step  in  designing  an  FIR  filter  is  to  build  the  vector  of  frequency  response  samples. 
For  a  low-pass  filter,  the  elements  which  belong  to  the  passband  of  the  ideal  fi'equency 
response  will  be  set  to  ‘1’  and  the  elements  which  belong  to  the  stopband  of  the  ideal 
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frequency  response  will  be  set  to  ‘O’.  The  number  of  elements  set  to  ‘1’  or  ‘0’  will  depend  on 
the  normalized  cutoff  frequency  and  on  whether  the  number  of  coefficients  is  odd  or  even. 
The  number  of  elements  to  set  to  ‘1  ’  in  the  frequency  response  vector  of  a  low-pass  filter  can 
be  computed  [5]  using  the  program  code  shown  in  Figure  28. 


if (num_taps%2 )  //odd  number  of  coefficients  case 

numsamples  =  ceil (num_taps  *  Wd/(2*PI)  -  0.293); 

else 

numsamples  =  ceil  (num__taps  *  Wd/{2*PI)  -  0.207); 


Figure  28  Computation  of  the  number  of  samples  to  include  in  the  passband  of  an 
ideal  low-pass  filter 

The  other  elements  of  the  vector  will  be  set  to  ‘0’  to  represent  the  stopband.  The  creation  and 
initialization  of  this  vector  is  equivalent  to  sampling  an  ideal  frequency  response  because 
each  element  of  this  vector  may  be  viewed  as  a  sample.  The  program  code  needed  to  build 
this  vector  is  shown  in  Figure  29.  Thus,  an  example  of  a  frequency  response  vector  could  be 
H[n\  =  [1,1,1,1,1,0,0, 0,0,0,]  for  a  low-pass  filter  with  10  coefficients. 


//sampling  of  the  ideal  frequency  response 

for{int  i=0;i<numsamples;i++) 

{ 

H[i]  =  1.0; 

} 

for ( int  i“numsamples ; i<num_taps ; i++ ) 

{ 

H[i]  =0.0; 

} 


Figure  29  Sampling  of  an  ideal  low-pass  filter  frequency  response 


4.3.2.1.2  Implementation  of  the  IDFT 

The  IDFT  is  performed  on  the  frequency  domain  samples  of  the  ideal  low-pass  filter 
generated  in  the  previous  section.  The  appropriate  expression  of  the  IDFT  is  used  depending 
on  whether  the  number  of  coefficients  is  odd  or  even.  The  algorithm  is  listed  in  Figure  30. 
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Float  temp  =0; 

Float  mid_pt  =  (num_taps-l . 0) /2 . 0; 

Float  x; 

if {num_taps%2)  //  N  odd 

{ 

for(int  n=0;n<num_taps;n++) 

{ 

temp  =  H [0] ; 

X  =  2  *  PI  *  (n  -  mid_pt) /num_taps; 
for (int  k=l; k< ( (num_taps-l) /2) ; k++) 
temp+= (2 . 0*cos (x*k) ) * (H [k] )  ; 
h[n]  =  temp/num_taps; 

fprintf  (outputfile,  ''%f \n",  h  [n] ) ;  / /printing  the 

coefficients  in  an  output  file 

} 

) 

else  //  N  even 

{ 

for (int  n=0;n<num_taps;n++) 

{ 

temp  =  H [0] ; 

X  =  2  *  PI  *  (n  -  mid_pt) /num_taps; 

for  (int  k=l;  k<=  (n\am_taps/2)  -  l;k++) 
temp+= (2 . 0*cos (x*k) ) * (H [k] )  ; 
h[n]  =  temp/num_taps; 
fprintf (outputfile, "%f \n", h [n] ) ; 

} 

) 


Figure  30  Inverse  discrete  Fourier  transform  algorithm  performed  on  the  ideal  low- 
pass  filter  frequency  response  samples 


4.3.2.1.3  Design  by  windowing 

The  only  difference  between  the  “Design  by  Windowing”  method  and  the  “Frequency 
Sampling  Design”  method  is  that,  for  the  design  by  windowing  method,  the  result  of  the 
IDFT  will  be  multiplied  by  a  vector  containing  the  samples  of  a  Hamming  window.  The 
equation  to  compute  the  samples  of  a  Hamming  window  is  given  in  Table  6  in  Section  4.2.3. 
The  computed  samples  of  the  Hamming  window  will  be  stored  in  a  vector  of  N  elements, 
where  N  is  the  number  of  filter  coefficients.  Then  the  coefficients  obtained  fi-om  the 
computation  of  the  IDFT,  and  stored  in  h[n],  will  be  multiplied  by  this  vector,  as  shown  is 
Figure  31. 
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for  (int  k=l;  k<  (  (nuiti__taps“l)  /2) ;  k++) 
temp+=  (2 . 0*cos  (x*k)  )  *  (H  [k]  )  ; 

//multiplying  the  time  domain  coefficients  by  the  Hamming  window 
h[n]  =  temp/num_taps  *  (0.54  -  0. 46*cos (2*PI*n/num_taps) ) ; 
//printing  the  coefficients  in  an  output  file 
fprintf (outputfile, "%f \n”, h [n] ) ; 


Figure  31  “Design  by  Windowing”  algorithm 

4.3.2.2  Parks-McClellan  method  implementation 

The  implementation  of  the  Parks-McClellan  FIR  coefficients  computation  module  is  included 
in  the  package  digital  Jitters  under  the  directory  named  remezex.  The  implementation  is 
separated  into  two  files:  remezex. c  and  remez.c.  ITie  library  remez.h  is  required  by  these  two 
files  and  contains  the  appropriate  constants  and  functions  used  by  them.  Hie  file  remez.c  and 
remez.h  are  external  files  that  were  first  created  by  Jake  Janovetz  and  can  be  used  under  the 
GNU  General  Public  license  restrictions.  The  source  code  implemented  by  Jake  Janovetz 
provides  a  very  good  implementation  of  the  Remez  Exchange  algorithm  and  forms  the  basis 
of  the  implementation  of  the  reconfigurable  module.  The  GNU  General  Public  License 
allows  anyone  to  use  and  modify  a  file  or  a  portion  of  code  placed  under  the  terms  of  this 
license.  TTie  GNU  General  Public  License  has  been  included  in  the  digital  Jitters  package. 
To  find  out  more  about  the  terms  and  agreements  of  this  license,  the  reader  is  referred  to 
http ://www. gnu.org/licenses/ gpl.html  on  GNU’s  web  site. 

The  task  of  the  file  remezex.  c  is  to  acquire  the  filter  design  parameters  entered  by  the  user, 
initialize  the  variables  used  by  the  module  and  make  the  flmction  calls  to  compute  the  filter 
coefficients.  It  is  noted  that  remezex.c  will  normalize  the  input  digital  frequency  over  the 
interval  [0,0.5],  as  is  required  for  the  implementation  of  the  Parks-McClellan  method  [3]. 
Then,  the  file  remezex.  c  fills  three  arrays  depending  on  the  design  parameters  entered. 

The  array  desired[]  contains  the  magnitude  of  the  sampled  ideal  frequency  response.  For  a 
low-pass  filter,  there  are  two  elements  in  the  array  desired[J  to  represent  the  magnitude  of  the 
fi’equency  response  in  the  passband  and  in  the  stopband  respectively.  In  this  implementation, 
the  first  element  is  initialized  to  ‘1’  to  represent  the  magnitude  of  the  frequency  response  in 
the  passband.  The  second  element  is  initialized  to  ‘0’  to  represent  the  magnitude  of  the 
response  in  the  stopband. 
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The  array  weights[]  contains  floating  point  values  representing  the  weights  given  to  the 
passband(s)  and  stopband(s).  These  weights  will  determine  the  importance  to  give  to  the 
corresponding  frequency  band  in  the  computation  of  the  algorithm.  A  large  value  given  to  a 
weight  corresponding  to  a  specific  band  will  put  the  emphasis  of  the  algorithm  on  producing 
a  better  approximation  of  the  ideal  frequency  response  in  this  band  compared  to  the 
approximation  made  in  the  other  bands  [1].  For  a  low-pass  filter,  the  weights []  array  contains 
two  elements.  The  first  one  specifies  Ae  weighting  assigned  to  the  passband  and  second 
specifies  the  weighting  assigned  to  the  stopband.  For  this  implementation,  each  frequency 
band  was  given  the  same  weighting  and  thus,  both  elements  are  set  to  ‘  1  ’ . 

The  array  bands[]  contains  the  edge  frequencies  which  delimit  the  passband  and  the 
stopband  of  the  low  pass  filter’s  frequency  response.  This  array  contains  four  normalized 
frequency  values  to  delimit  the  passband  and  the  stopband  and  should  be  in  the  form ,  [0,  Op, 
(Oa,  0.5].  In  this  implementation,  the  user  enters  the  desired  normalized  cutoff  frequency,  Wd, 
and  a  value  of  Q)p=Wj- 0.025  and  (o^=Wj+  0.025  is  computed  for  the  array  bands[]. 

The  length  of  the  transition  band  between  the  passband  and  the  stopband,  for  this  module, 
has  been  arbitrarily  chosen  to  be  0.05  Hz/Hz. 

Figure  33  shows  the  contents  of  remezex.c.  After,  filling  the  three  required  arrays,  remezex.c 
makes  a  fimction  call  to  void  remezQ,  the  ftmction  which  computes  the  filter  coefficients 
using  the  Remez  Exchange  algorithm,  void  remezQ  is  implemented  in  the  file  remezx. 


desired[0] 

=  1; 

desired [1] 

=  0; 

weights [0] 

=  1; 

weights [1] 

=  1; 

bands [ 0 ]  = 

0; 

bands [ 1 ]  = 

Wd-0.025; 

bands [2]  = 

(Wd  +  0.025) ; 

bands [3]  = 

0.5; 

remez ( &h [0] ,nuin_taps,  2,  bands,  desired,  weights, BANDPASS) ; 

Figure  32  Content  of  the  file  remezex.c 


4.3.2.3  Gaussian  filter  implementation 

The  implementation  of  the  Gaussian  filter  coefficients  computation  module  is  based  on  the 
“Frequency  Sampling  Design”  method.  Instead  of  computing  the  samples  of  the  frequency 
response  of  an  ideal  low-pass  filter,  samples  for  a  Gaussian  distribution  are  calculated  based 
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on  Equation  (98).  Figure  32  shows  the  portion  of  code  that  is  used  to  fill  in  the  array  of 
frequency  response  samples,  H[i]. 


//computing  the  sampling  step 

step  =  (float) (Fs/num_taps) ; 


//sampling  of  the  gaussian  curve 
for(int  i«=0;i<num_taps;  i++) 

( 

H[i]  =  exp( (-1.0/4. 0)*(2.0*PI*i*step)* 

(2 . 0*PI*i*step) / (K*K*B*B) ) ; 


} 


Figure  33  Sampling  algorithm  of  a  Gaussian  distribution 

The  next  step  of  the  implementation  is  to  compute  the  IDFT  of  the  sampled  curve  to  get  the 
filter  coefficients  as  was  done  for  the  “Frequency  Sampling  Design”  method  implementation 
and  the  “Design  by  Windowing”  method  implementation. 


4.3.3  FIR  signal  processing  modules 

The  FIR  signal  processing  module  has  been  implemented  based  on  the  convolution  equation 
introduced  in  Section  4.1.  This  implementation  is  based  on  the  Direct  Form  structure 
presented  in  Section  4.1.1.  The  signal  processing  modules  can  be  found  in  the  package  of 
filter  modules  under  the  name  offirdf. 

Figure  34  shows  the  implementation  of  the  algorithm  of  the  FIR  signal  processing  module.  A 
buffer  is  used  to  store  the  input  samples  required  for  the  convolution.  First,  the  previous  input 
samples  are  shifted  down  in  the  buffer.  The  current  input  sample  is  then  read  from  an  input 
text  file  and  stored  at  the  beginning  of  the  buffer.  Then,  a  convolution  sum,  as  described  in 
Equation  (4),  is  performed  to  produce  an  ouq)ut  sample  which  is  stored  in  an  array.  The 
output  array  is  printed  to  a  text  file  once  all  the  input  samples  have  been  processed.  This 
implementation  is  able  to  operate  in  real  time  as  each  input  sample  is  read  and  processed  to 
directly  produce  an  ou^ut  sample. 
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while (fscanf (input file, "%f ”, &sample) !=EOF) 

//reading  input  sample  from  the  input  file 

for ( j=numcoef-l; j>0; j — ) 

{ 

if{indx-j  <  0)  continue; 

inputBuf fer [ j ] =inputBuf fer [ j-1]  ; 

} 

inputBuf fer [0] =sample; 

//convolution  operation  to  produce  output  sample 
output=0; 

for(i  =  0;  i<  numcoef;I++) 

{ 

if{indx-i  <  0)  break; 
output+=coeff [i] *inputBuf fer [i] ; 

} 

fprintf (outputfile, ”%f\n”, output ) ; 
indx++; 

} _ 


Figure  34  FIR  signal  processing  module  algorithm 


4.3.4  Use  of  the  filter  and  signal  processing  modules 

The  following  is  an  example  of  how  the  modules  are  used  to  generate  coefficients  for  a  low- 
pass  Butterworth  filter.  The  coefficients  computation  modules  can  be  executed  from  the 
prompt  of  a  console  as  shown  in  Figure  35.  The  module  butter  takes  four  arguments  as  input 
design  specifications,  the  filter’s  order,  the  normalized  cutoff  frequency,  the  sampling 
frequency  and  the  number  associated  with  the  method  used  to  compute  the  discrete  design. 
The  filter  order  is  entered  as  an  integer  value.  The  normalized  cutoff  frequency  is  entered  as 
a  floating  point  number  between  0.0  and  1.0  for  the  Butterworth  filter.  The  sampling 
frequency  is  a  floating  point  number  and  is  normalized  to  2.0  in  this  example.  The  units  of 
the  normalized  cutoff  frequency  and  the  sampling  frequency  are  Hz/Hz.  It  is  noted  that  the 
normalized  cutoff  and  sampling  frequencies  are  relative.  As  a  result,  the  filter  module  allows 
for  scalability  of  the  filter  design  and  generates  the  same  coefficients  for  filters  with  the  same 
ratio  of  cutoff  frequency  to  sampling  frequency.  The  last  parameter  specifies  the  method 
used  to  compute  the  discrete  design.  A  ‘1’  is  entered  to  select  the  Impulse  Invariance 
method,  and  a  ‘2’  is  entered  to  select  the  Bilinear  transformation. 
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After  the  computation  of  the  coefficients  is  complete,  the  module  generates  an  ou^ut  file, 
“butterout.txf’  with  the  coefficients  of  the  filter.  Figure  36  shows  an  example  of  the 
coefficients  computation  module  output  file.  The  order  of  the  filter  computed  is  printed  on 
the  first  line.  Then,  on  the  second  line,  the  digital  static  gain  is  printed  and  finally,  the 
coefficients  of  the  filter  are  printed  in  four  columns.  The  two  left  most  columns  contain  the 
real  and  imaginary  parts  of  the  numerator  coefficients  (6*).  The  two  right  most  columns 
contain  the  real  and  imaginary  parts  of  denominator  coefficients  (a*). 


1  ^  bulterout.tKt  -  Notepad 

JPllxt 

Ij  File  Edit  Format  Help 

Mi 

5.000000  - 

0.052786  - 

0.052786  0.000000  1.000000  0.000000 
0.263932  0.000000  0.000000  0.000000 
0.527864  0.000000  0.633437  0.000000 
0.527864  0.000000  0.000000  0.000000 
0.263932  0.000000  0.055728  0.000000 
0.052786  0.000000  0.000000  0.000000 

I  T^fl>nrttnin!>+r 


Order 


Denominator  Coefficients 
Numerator  Coefficients 


Figure  36  Example  of  an  output  file  generated  by  the  coefficients  computation 

modules. 
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The  signal  processing  modules  can  be  executed  from  a  console  prompt  as  in  the  case  of  the 
coefficients  computation  modules.  Figure  37  shows  an  example  of  how  to  execute  the  signal 
processing  module  for  an  HR  filter.  The  first  argument  of  these  modules  is  the  file  name  of 
the  input  samples.  The  second  argument  is  the  name  of  the  text  file  which  contains  the 
coefficients  of  the  implemented  filter  design.  The  coefficients  file  is  the  output  file  obtained 
from  the  coefficients  computation  modules  (i.e.'^butterouf).  The  user’s  guide  contains 
information  on  the  format  of  the  input  data  file. 


Command  Prompt 


□ax 


Microsoft  U in do us  2G6C  LUers ion  C-CG.2i9B] 


(C>  i9CS-2G<jC  [ 'ici'ost;  i' ;;  . 

C:\>iil^  ScJToiecf  i  le  cooffile,  ^ 


Samples  flic 
-  C'ociricients  file 


n 


Figure  37  Command  line  for  executing  an  HR  signal  processing  module. 


The  output  file  generated  by  the  signal  processing  module  contains  the  output  of  the  filter 
used  to  process  the  input  samples  file.  One  output  sample  is  printed  on  each  line  of  the  file  as 
shown  in  Figure  38.  The  first  column  contains  the  real  part  of  the  output  and  the  second 
column  contains  the  imaginary  part  of  the  output.  In  the  example  shown  in  Figure  38,  the 
values  of  the  imaginary  part  of  the  results  are  zero  as  a  real  signal  was  processed  rather  than  a 
complex  signal. 


lIRoutput.txt  -  Notepad 


File  Edik  Format  Help 


l£l|x| 


0.178370 

0.000000 

-0.160113 

0.000000 

-0,184  593 

0.000000 

0.053849 

0.000000 

0.137787 

0.000000 

-0.017965 

0.000000 

-0.098050 

0.000000“ 

0.005989 

0.000000 

0. 069173 

0.000000 

-0.001996 

0. 000000 

-0.048722 

0.000000 

0.000665 

0.000000 

Jj 

-  Real  Part 
Imaginary  Part 


Figure  38  Example  of  an  output  file  generated  by  a  signal  processing  module 
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5  DSP  IMPLEMENTATION 


The  implementations  of  the  filter  modules  described  in  the  previous  sections  have  been 
targeted  for  a  general  purpose  processor.  In  order  to  make  use  of  these  digital  filter  modules 
in  the  ROBR  and  take  advantage  of  their  reconfigurability  in  “real-time”,  a  DSP 
implementation  of  the  modules  is  required.  For  the  DSP  implementation,  a  DSP  board 
manufactured  by  Spectrum  Signal  Processing,  Inc.,  was  used  The  DAYTONA  DSP  board  is 
a  PCI  dual  processor  board  and  contains  two  TMS320C6201  fixed  point  DSP  chips.  Only 
one  processor  is  used  for  the  filter  module  implementation.  Code  development  for  the  DSP 
board  is  done  in  the  ANSI  C  programming  language.  A  compiler,  code  generator,  and  linker 
are  provided  with  the  DSP  board. 

Two  FIR  coefficient  computation  modules  and  one  FIR  signal  processing  module  have  been 
adapted  for  the  DAYTONA  DSP  board.  In  all  three  cases,  since  the  original  module  was 
written  in  C,  only  minor  changes  were  required  to  yield  compatibility  between  the  filter 
modules  and  the  DSP  architecture.  The  coefficients  computation  modules  implemented  are 
the  “Design  by  Windowing”  method  and  the  “Parks-McClellan”  method.  The  FIR  signal 
processing  module  implements  the  Direct  Form  structure  shown  previously  in  Section  4.1.1. 

The  DSP  board  is  hosted  in  a  PC  with  the  WIN  NT  operating  system.  The  host  handles  the 
initialization,  handshaking  and  downloading  of  processor  code  through  the  PCI  bus.  As 
such,  the  integration  of  the  filter  and  signal  processing  modules  requires  two  different  files:  a 
host  program  that  controls  the  DSP  and  provides  the  user  interface,  and  the  DSP  file  which 
implements  the  DSP  program  is  responsible  for  processing  the  data. 


5.1  Exchanging  data  between  the  Daytona  and  the  host  station 


The  DSP  processor  has  both  internal  and  external  memory  spaces  available  including  internal 
program  and  data  RAM,  external  SSRAM,  external  SDRAM,  and  dual  port  RAM.  For  this 
project,  the  memory  used  for  the  data  exchange  is  the  SDRAM  of  processor  0  on  the  Daytona 
DSP  board  [9].  The  SDRAM  block  of  the  processor  goes  from  adless  0x0200000  to  address 
0x02FFFFFF.  The  other  memory  spaces  available  on  the  DSP  board  are  listed  in  Table  8. 


Description 

Size 

Internal  or  External 

Program  RAM 

64kB 

Internal 

SSRAM 

16MB 

External 

I/O,  boot 

4MB 

External 

SDRAM 

16MB 

External 

Processor  Expansion 

Module  (PEM) 

64kB 

External 

Internal  Registers 

256kB 

Internal 

Table  8  Memory  configuration  of  the  TMS320C6201 
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5.2  Static  and  determined  length  memory  allocation 


The  static  allocation  of  memory  is  done  on  the  SDRAM  block  of  the  processor.  Figure  39 
shows  the  memory  allocation  for  a  typical  filter  module.  Four  32-bit  variables  are  used  as 
flags  for  each  module.  They  are  Flag_Ready,  Flag  Data,  Flag_Done,  and  OK  Memory. 
The  addresses  for  the  memory  space  are  assigned  in  bytes  so  that  consecutive  32-bit  words 
are  addressed  4  bytes  apart.  The  use  of  flags  makes  die  synchronization  between  the  host 
program  and  the  DSP  possible.  The  flag  Flag_Ready  is  set  by  the  DSP  to  tell  the  host 
program  that  the  DSP  is  ready  to  compute  data.  The  flag  Flag  Data  is  set  by  the  host  to 
signal  to  the  DSP  that  data  is  ready  to  be  transferred  from  the  host.  The  flag  Flag_Done  is  set 
by  the  DSP  to  tell  the  host  program  that  the  DSP  has  finished  its  computations.  The  flag 
OK_Memory  is  used  by  the  DSP  to  tell  the  host  program  that  the  memory  needed  for  the 
computation  has  been  successfully  allocated.  As  shown  in  Figure  39  memory  has  been 
allocated  for  Num  Input  and  NumjCoeff  which  are  integer  values  that  are  used  to  store  the 
number  of  coefficients  to  be  computed  and  the  number  of  input  data  samples.  Arrays  are 
allocated  at  addresses  pointed  to  by  Coeff  and  Data.  The  array  Coeff  contains  the  coefficients 
computed  by  the  DSP  that  are  transferred  back  to  the  host  program.  The  array  Data  is  used  to 
hold  the  input  signal  samples  to  be  transferred  from  the  host  to  the  DSP.  After  the  DSP  has 
performed  the  signal  processing  computation,  the  same  area  is  used  to  store  the  computed 
output  signal.  The  host  program  will  then  be  able  to  retrieve  the  ou^ut  signal  from  the  array 
Data.  Figure  40  shows  a  flow  diagram  of  the  handshaking  between  Ae  host  and  the  DSP. 


#define  Flag_Ready  (UINT32* ) (0x02000004 ) 
#define  Flag_Data  (UINT32*) (0x02000008) 
#define  Flag_Done  (UINT32*) (0x0200000c) 
#define  Num_Input  (UINT32*) (0x02000010) 
♦define  Nuin_Coeff  (UINT32*)  (0x02000014) 
♦define  OK_Memory  (UINT32*) (0x02000018) 
♦define  Coeff  (float*) (0x0200001c) 
♦define  Data  (float*) (0x0200090c) 


Figure  39  Static  allocation  of  the  variables  in  the  SDRAM  memory 


Host 


Figure  40  Flow  diagram  of  the  implementation  of  the  filter  modules 


5.3  Dynamic  memory  allocation 

The  function  call  mallocQ  is  a  service  provided  by  the  run-time  support  library  included  in 
the  DSP  compiler.  mallocO  extensively  used  in  the  reconfigurable  modules  because  the  size 
of  the  arrays  needed  to  store  the  data  from  the  user  is  not  known  in  advance.  Memory  is 
dynamically  allocated  from  a  memory  space  defined  in  the  “sysmem”  memory  section  of  the 
DSP.  The  sysmem  memory  section  is  created  and  allocated  prior  to  the  compilation  and  the 
linking  of  the  program  in  Ae  link  command  file.  In  Figure  41,  a  portion  of  the  memory  map 
generated  by  the  linker  shows  that  the  “sysmem”  memory  section,  which  is  0x3000  bytes  (or 
12kB)  in  size,  is  allocated  to  begin  at  address  0x80002018. 
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TMS320C6X  COFF  Linker  Version  2.00 


MEMORY  CONFIGURATION 

name  origin 

length 

used 

attributes  fill 

IVECS 

00000000 

000000400 

00000200 

RWIX 

I  PROG 

00000400 

OOOOOfcOO 

OOOOlbSO 

RWIX 

SSRAM 

00400000 

000400000 

00000000 

RWIX 

MPRAM 

01400000 

000200000 

00000000 

RWIX 

DL3 

01600000 

OOOOcOOOO 

00000000 

RWIX 

I  REG 

01800000 

000800000 

00000000 

RWIX 

SDRAM 

02000000 

001000000 

00000000 

RWIX 

PEM 

03000000 

001000000 

00000000 

RWIX 

IVARS 

80000000 

OOOOcOOOO 

000050e8 

RWIX 

I  DATA 

800C0000 

000004000 

00000000 

RWIX 

SECTION 

ALLOCATION  MAP 

output 

attributes/ 

section 

page 

origin 

length 

input  sections 

.vectors 

0 

00000000 

00000200 

00000000 

00000200 

isfp6201.o6x  (.vectors) 

.text 

0 

00000400 

OOOOlbOO 

00000400 

000009C0 

rts6201.1ib  ;  memory. obj  (.text) 

.stack 

0 

80000000 

00002000 

UNINITIALIZED 

80000000 

00000000 
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Figure  41  Sysmem  memory  section  allocation  in  the  DSP  memory  map 


5.4  Host  program 


The  major  function  of  the  host  program  is  to  control  the  DSP  board  and  act  as  the  user 
interface  to  the  DSP  board.  The  host  program  is  responsible  for  reading  the  design 
specifications  of  the  filter  module  and  passing  them  to  the  DSP.  The  Daytona  Windows  NT 
Host  Application  Library  (ALIB_HOST)  provides  several  high-level  functions  that  allow  the 
user  to  control  the  operations  of  the  Daytona  from  a  Windows  NT  host. 

5.4.1  Host  software  functions 

The  library  ALIB_HOST  implements  several  host  functions  to  control  the  Daytona  board. 
The  following  functions  have  been  used  in  the  filter  modules’  implementation. 


FT_ControlO 

FT_ErrorMessageO 

FT_GetHandle() 

FT_ReadO 

FT_SystemClose() 

FT_SystemLoadO 

FTSystemOpenO 

FT_Write() 


:  to  reset  the  board 
:  to  catch  error  message 
:  to  get  handle  to  Daytona  system  resources 
:  to  read  fix)m  a  system  resource  or  host  buffer 
:  to  close  the  DSP  board 
:  to  load  the  DSP  code  into  the  system 
:  to  open  the  system 
:  to  write  to  the  system  resources 


5.5  DSP  program 

The  DSP  program  contains  all  the  data  processing  instructions.  The  algorithm  of  the  original 
modules  has  not  been  modified  but  the  variables  declared  were  changed  into  pointers  in  most 
cases  to  be  integrated  for  use  on  the  DSP. 

5.5.1  DSP  software  functions 


As  in  the  host  program,  library  functions  are  provided  for  initialization,  interrupts,  and  DMA 
transfers.  In  this  project,  the  only  library  function  called  in  the  DSP  code  is  C6x_OpenC6xO. 
This  function  initializes  the  C6x  processor,  sets  the  wait  states  for  external  memory  and 
configures  the  page  register.  The  page  register  contains  the  addresses  of  the  memory  spaces 
that  are  available  for  each  processor.  The  memory  spaces  are  accessed  through  the  PCI  bus 
[11]. 
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6  Results  and  Verification 
6.1  Methodology 

The  implementation  of  the  filter  modules  has  been  verified  with  the  Signal  Processing 
toolbox  from  MATLAB.  Table  9  and  Table  10  show  which  corresponding  functions  from 
MATLAB  have  been  used  to  verify  the  coefficients  generated  by  the  reconfigurable  filter 
modules. 


HR  coefficient  computation  modul 

es 

Filter  types 

Reconfigurable  Filter 
modules 

Function  from  MATLAB’s 
toolbox  signal 

Butterworth 

butter 

Butter 

Chebyshev 

chebl 

chebyl,  cheblap 

Inverse  Chebyshev 

cheb2 

cheby2,  che2ap 

Elliptical 

ellip 

ellip,  ellipap 

Bessel 

bessel 

- 

Table  9  Corresponding  MATLAB  functions  for  the  HR  modules  verification 


FIR  coefficient  computation  modules 

Design  method  used  and 
Filter  types 

Reconfigurable  Filter 
modules 

Function  from  MATLAB’s 
toolbox  signal 

Frequency  sampling 

fi'eqsampling 

- 

Windowing 

hamming 

firl  with  hamming  window 

Parks-McClellan  (optimal 
equiripple  filter) 

remezex 

Remez 

Gaussian  filter 

gauss 

- 

Digital  Integrator 

dintegrator 

- 

Table  10  Corresponding  MATLAB  functions  for  the  HR  modules  verification 


The  signal  processing  modules  for  the  HR  and  FIR  filters  were  compared  with  the 
corresponding  yi/rer  function  fi'om  the  Signal  Processing  toolbox.  Coefficients  generated  by 
the  coefficient  computation  modules  were  provided  to  the  appropriate  signal  processing 
module  to  compute  the  response  to  an  impulse.  Subsequently,  an  FFT  is  performed  using 
MATLAB ’s  function  on  the  output  of  the  signal  processing  module  which  gives  the 
frequency  response.  Samples  of  the  impulse  function  are  used  as  the  input  to  the  signal 
processing  module.  The  samples  are  generated  using  MATLAB  and  are  composed  of  49 
“zeroes”  followed  by  a  “one”,  followed  by  49  “zeroes”  in  floating  point  representation. 
Figure  42  is  a  plot  of  the  impulse  signal  generated  with  MATLAB. 
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1 


Test  signal 


Figure  42  Impulse  signal  generated  with  Matlab 

Thus,  the  verification  process  was  carried  out  as  follows: 

1-  Computation  of  the  coefficients  with  the  coefficients  computation  modules 

2-  Computation  of  the  coefficients  with  MATLAB ’s  toolbox 

3-  Comparison  of  the  coefficients  generated  by  the  reconfigurable  modules  with  the 
coefficients  generated  by  MATLAB 

4-  Processing  of  an  impulse  with  the  signal  processing  modules  using  the  coefficients 
computed  by  the  coefficients  computation  modules 

5-  Processing  of  an  impulse  with  the  MATLAB’s  function  using  the  coefficients 
computed  by  MATLAB 

6-  Performing  an  FFT  on  both  output  generated  by  the  filter  modules  and  by  MATLAB 
functions. 

7-  Comparison  of  the  frequency  response  yielded  by  the  filter  modules  and  by 
MATLAB. 


6.2  Results 

6.2.1  HR  filter  module  verification 

Appendix  A  shows  the  plots  comparing  the  frequency  response  of  the  Butterworth, 
Chebyshev,  inverse  Chebyshev,  and  elliptical  filter  modules  with  those  generated  by 
MATLAB.  The  frequency  response  for  the  bessel  filter  module  is  also  given.  However,  there 
is  no  implementation  of  a  digital  bessel  filter  in  MATLAB.  In  all  cases,  the  responses  show 
a  steeper  rolloff  for  the  higher  order  filters.  Both  coefficients  computed  by  the  filter  modules 
and  by  MATLAB’s  functions  are  identical  for  filter  orders  up  to  15.  In  cases  of  more  than  16 
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coefificients,  the  accuracy  of  the  frequency  response  yielded  by  the  coefiBcients  starts  to 
decrease.  This  limitation  may  be  caused  by  the  rounding  off  of  floating  point  number 
variables  with  the  large  number  of  multiplication  and  addition  operations  used  to  compute 
the  filter  coefficients. 

6.2.2  FIR  filter  module  verification 

For  FIR  filters,  a  higher  order  filter  corresponds  to  a  higher  number  of  coefficients.  From 
Appendix  A,  the  three  FIR  filter  design  methods  show  that  as  the  number  of  coefficients 
increases,  the  transition  between  the  passband  and  stopband  is  much  steeper,  as  expected.  To 
illustrate  this  characteristic,  an  error  curve  is  plotted  showing  the  difference  between  the 
frequency  response  of  an  ideal  low-pass  filter  and  the  frequency  response  of  the 
reconfigurable  filter  for  each  filter  design  method.  An  example  of  an  error  curve  for  the 
“Frequency  Sampling  Design”  method  with  N=10  and  N=30  coefficients  is  shown  in  Figure 
43. 
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Figure  43  Error  curve  for  reconfigurable  filter  using  the  “Frequency  Sampling 
Design”  method,  (a)  for  N=10  coefficients,  (b)  for  N=30  coefficients. 
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The  spike  in  the  error  curve  corresponds  to  the  transition  between  the  passband  and  stopband. 
It  is  noted  that  for  N=30,  the  spike  is  more  compressed  along  the  x-axis  than  for  the  N=10 
case,  indicating  a  faster  transition  between  the  passband  and  stopband. 

As  with  the  IIR  filter  modules,  the  frequency  responses  of  the  reconfigurable  FIR  filters  are 
also  presented  in  Appendix  A.  The  power  spectrum  for  each  of  these  filters  is  also  shown. 

In  all  cases,  the  FIR  filters  are  designed  to  have  a  normalized  cutoff  frequency  of  0.5. 

The  FIR  coefficients  computation  modules  that  use  the  “Design  by  Windowing”  method  and 
the  “Parks-McClellan”  method  have  been  verified  with  MATLAB’s  respective  functions  and 
are  shown  in  Appendix  A.  It  is  shown  that  the  reconfigurable  filter  modules  yield  the  same 
general  shape  of  frequency  responses  as  their  respective  counterpart  in  MATLAB.  The  FIR 
coefficients  computation  modules  for  the  “Frequency  Sampling  Design”  method,  the 
Gaussian  filter  coefficients  computation  module  and  the  digital  integrator  have  not  been 
verified  with  MATLAB’s  toolbox  because  the  corresponding  functions  in  MATLAB  were 
not  available. 

As  shown  in  the  error  curves,  the  “Frequency  Sampling  Design”  method  yields  frequency 
responses  for  which  the  transition  between  passband  and  stopband  becomes  steeper  as  the 
filter  order  increases.  A  consequence  of  the  steeper  transition  is  an  overshoot  of  the 
frequency  response  at  the  start  of  the  transition  which  measured  approximately  1 1%  of  the 
passband  magnitude.  The  attenuation  in  the  stopband  is  shown  to  start  at  — 15dB  and 
gradually  rolls  off  to  approximately  -30dB. 

The  “Design  by  Windowing”  method  yields  better  frequency  responses  with  less  ripple  in  the 
stopband  and  in  the  passband,  and  sm^ler  overshoot  at  the  transition.  Furthermore,  the 
results  indicated  significantly  better  attenuation  in  the  stopband  at  around  -50dB.  However, 
when  observing  the  error  curves  in  Table  AlO  of  Appendix  A,  a  wider  spike  suggests  a  more 
gradual  transition  between  the  passband  and  stopband. 

The  frequency  responses  of  the  “Parks-McClellan”  filters  exhibited  ripples  in  both  the 
passband  and  stopband  as  did  the  “Frequency  Sampling  Design”  method.  A  difference  in  the 
amplitude  of  the  frequency  responses  between  the  reconfigurable  filter  module  and  the 
MATLAB  functions  can  be  attributed  to  different  values  of  the  passband  and  stopband  ripple. 
With  a  greater  stopband  ripple  allowed  in  the  MATLAB  case,  less  attenuation  is  noted  in  the 
power  spectrum  when  compared  with  the  reconfigurable  filter  modules.  The  error  curves  for 
the  Parks-McClellan  method  demonstrated  transitions  between  passband  and  stopband  that 
were  comparable  to  the  “Frequency  Sampling  Design”  method.  While  the  Parks-McClellan 
method  provides  flexibility  in  setting  the  passband  and  stopband  ripple,  it  is  much  more 
computationally  intensive.  This  may  have  sigmficant  implications  when  choosing  a  filter 
design  method  for  real-time  processing  requirements. 

Frequency  response  and  power  spectrum  curves  are  also  plotted  for  a  Gaussian  filter  with 
BT=0.2  for  N=10,20,  and  30  coefficients.  As  mentioned  in  Section  4.3.2.2,  the  “Frequency 
Sampling  Design”  method  was  used  to  implement  the  reconfigurable  Gaussian  filter  module. 


62 


The  results  showed  that  while  the  main  lobe  did  not  change  significantly  as  the  number  of 
coefficients  increased,  more  atteniiation  in  the  stopband  was  observed  as  N  increased. 

Results  for  the  digital  integrator  module  are  also  shown  in  Appendix  A.  The  output  curve 
shows  the  integration  of  a  bipolar  input  digital  bit  stream. 
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7  Summary 


This  report  describes  a  project  to  develop  reconfigurable  HR  and  FIR  filter  modules  for  use 
in  the  ROBR  project.  The  theory  needed  to  understand  digital  signals,  LTI  systems  and  filters 
was  introduced  in  Section  2.  llie  concept  of  HR  and  FIR  filters  was  presented  in  Section  3 
and  Section  4.  Both  the  HR  and  FIR  implementations  were  represented  using  the  Direct  Form 
structure.  Equations  for  the  analog  frequency  response  of  various  types  of  HR  filters  were 
presented.  Examples  include  Butterworlii,  Chebyshev,  and  Elliptical  filters. 

The  design  methods  used  to  implement  HR  filters  produce  a  discrete  filter  design  fi'om 
analog  design  prototypes.  Two  methods  have  been  used  to  implement  HR  filters:  the  Bilinear 
transformation  and  the  Impulse  Invariance  method.  These  two  methods  perform  a 
transformation  on  an  analog  design  to  obtain  a  discrete  design. 

Methods  used  to  produce  FIR  filters  involve  sampling,  IDFT  computation  and  an 
optimization  algoritiim.  Three  design  methods  have  been  used  to  implement  FIR  filter 
modules:  the  “Frequency  Sampling  Design”  method,  the  “Design  by  Windowing”  method 
and  the  “Parks-McClellan”  me&od.  The  “Frequency  Sampling  Design”  method  involves  two 
steps  in  the  computation  of  the  filter  coefficients:  building  an  ideal  frequency  response 
vector,  and  computing  an  IDFT.  A  vector  is  a  discrete  sequence  of  elements.  Building  the 
vector  is  equivalent  to  sampling  the  ideal  frequency  response.  The  “Design  by  Windowing” 
method  involves  the  same  steps  as  the  “Frequency  Sampling  Design”  method  except  that  the 
result  of  the  IDFT  performed  on  the  desired  fi^quency  response  samples  are  subsequently 
multiplied  by  a  vector  containing  amplitude  samples  of  a  window  Action.  A  Hamming 
window  has  been  used  for  this  implementation.  The  “Parks-McClellan”  method  is  based  on 
the  Alternation  Theorem  firom  optimization  theory.  The  Remez  Exchange  algorithm  is  used 
to  find  the  optimal  set  of  extremal  fi'equencies.  The  goal  of  the  method  is  to  compute  the 
coefficients  for  the  best  approximation  of  a  desired  frequency  response.  The  Remez 
exchange  algorithm  is  a  set  of  conditional  statements  that,  when  applied,  produce  an  optimal 
fi'equency  response.  The  design  of  a  Gaussian  filter  using  the  “Frequency  Sampling  Design” 
me&od  is  also  presented. 

Reconfigurable  filter  modules  were  implemented  for  five  types  of  HR  filters  and  four  types 
of  FIR  filters.  The  implemented  types  of  HR  filters  are  the  Butterworth  filter,  the  Chebyshev 
filter,  the  Inverse  Chebyshev  filter,  the  Elliptical  filter  and  the  Bessel  filter.  Two  FIR  filter 
modules,  based  on  the  “Frequency  Sampling  Design”  method,  were  implemented  called 
freqsampling,  and  gauss,  the  latter  of  which  is  the  Gaussian  filter  module  implementation. 
An  FIR  filter  module  was  also  implemented  using  the  “Design  by  Windowing”  method  and 
is  called  hamming.  The  “Parks-McClellan”  FIR  implementation  is  called  remezex. 

Coefficients  generated  by  the  reconfigurable  filter  modules  were  used  to  process  an  input 
impulse  function.  The  resulting  firequency  response  was  compared  to  that  generated  by 
MATLAB  and  is  presented  in  Appendix  A.  The  digital  reconfigurable  HR  filter  modules 
were  found  to  yield  the  same  responses  as  the  corresponding  functions  in  MATLAB.  The 
FIR  digital  filter  modules  also  yielded  the  same  fi-equency  responses  as  MATLAB.  The 
Parks-McClellan  method  provides  the  flexibility  to  adjust  die  passband  and  stopband  ripple 
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while  yielding  a  comparatively  steep  transition.  However,  the  algorithm  is  much  more 
computationally  intensive.  A  Gaussian  filter  module  and  a  digital  integrator  module  were 
also  successfully  implemented.  The  Gaussian  filter  and  digital  integrator  modules  will 
facilitate  the  development  of  a  GMSK  modulator  for  the  ROBR. 

The  “Design  by  Windowing”  and  “Parks-McClellan”  filter  modules  were  successfully 
adapted  for  use  on  a  DSP  board.  As  well,  a  module  for  processing  a  signal  using  the  FIR 
filter  coefficients  generated  was  implemented  on  the  DSP.  The  DSP  board  used  was  the 
Daytona  Dual  c62  processor  board  from  Spectrum  Signal  Processing  Inc.  Issues  related  to 
dynamic  memory  allocation  and  handshaking  between  the  host  and  the  DSP  were  discussed. 
Further  work  is  required  to  completely  adapt  the  filter  modules  for  the  ROBR,  to  assess  the 
ability  to  reconfigure  the  modules  while  the  ROBR  is  operating,  and  to  resolve  any  time 
critical  issues  for  the  ROBR  where  more  computationally  intensive  algorithms  are  used. 
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Appendix  A 


A1.  Results 


A1.1.  (IR  filter  modules 


A1.1.1.  Butterworth  filter 


Table  A1  Results  for  the  Butterworth  filter  module 


A1 


A1.1.2.  Chebyshev  filter 


Table  A2  Results  for  the  Chebyshev  filter  module 
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A1.1.3.  Inverse  Chebyshev  filter 


Table  A3  Results  for  the  Inverse  Chebyshev  filter  module 
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A1.1.4.  Elliptical  filter 


Reconfigurable  filter  module 
firequency  responses 

Corresponding  MATLAB  fimction 
frequency  responses 
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Table  A4  Results  for  the  Elliptical  filter  module 
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A1.1.5.  Bessel  filter 


Reconfigurable  filter  module 
Frequency  responses 


Order 


Corresponding 
MATLAB  function 
frequency  responses 
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Table  A5  Results  for  the  Bessel  filter  module 


A1.2.  FIR  filtering  modules 

A1.2.1.  Frequency  sampling  design  method 
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Table  A6  Results  for  the  Frequency  Sampling  Design  filter  module 
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Table  A7  Frequency  and  power  spectrum  for  N=10  FIR  filter  using  the  Frequency 

Sampling  Design  filter  module 
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Table  A9  Frequency  and  power  spectrum  for  N=30  FIR  filter  using  the  Frequency 

Sampling  Design  filter  module 


A1.2.2.  Design  by  windowing  with  a  Hamming  window 


Filtering  modules 
coefficients 

Mean 

error 

Error  Curve 
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Reconfigurable  filter) 
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Table  A13  Frequency  and  power  spectrum  for  N=30  FIR  filter  using  the  Design  by 
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A1.2.3.  Parks-McClellan  method 


Filtering  modules 
coefficients 


Mean  Error  Curve 
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Table  A15  Frequency  and  power  spectrum  for  N=10  FIR  filter  using  the  Parks-McClellan 

filter  module 


A15 


Table  A16  Frequency  and  power  spectrum  for  N=21  FIR  filter  using  the  Parks-McClellan 

filter  module 


A16 


IHhI 


Table  A17  Frequency  and  power  spectrum  for  N=30  FIR  filter  using  the  Parks-McClellan 
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Steps  to  calculate  poles  and  zeros  of  the  elliptic  filter 


4 


The  Elliptical  filter  coefficients  computation  module  implements  the  following  steps  to 
calculate  the  poles  and  zeros  of  such  a  filter.  Poles  and  zeros  are  then  used  to  compute  the 
coefficients  of  the  filter  [4], 

(Op  =  passband  fi'equency 
fl),  =  stopband  frequency 
Ap  =  maximum  passband  loss  (dB) 

-  maximum  stopband  loss  (dB) 
k  =  selectivity  factor  =(Opl(o^ 


1.  k'  =Vl-it^ 


2.«.  =  2 


1-VF 

i+VF 


) 


3.  ^  =  ^0 +2^0* +15^0* 


c  ^  logl6D 

D.  n> — = - 

log(l/^) 


6.  A  = 


_1_  10°°'^^ 
2„“io0.o5^, 


-1 


7. 


0-0  = 


2q^'*  ^  (-1)"  ^ sinh[(2m  +  1)A] 

m»0 

l  +  2j(-l)'”^'"' cosh2wA 

OT*=1 


8.  IF  =1(1+ 


1  +  ^ 

V  ^  7 
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9.  n.  = 


29'"' 2  (-1)'"^”^'”"'^  sin 
«*0 

(2m  + 

n 

i+2j;(-ir^'"’cosh 

ImTtn 

n 

/'  = 


/ 

*  2 


For  n  odd 
For  n  even 


i  =  1,2,...,/ 


10.  F,=  + 


Q 


6o.= 


(o-oF,)^ 


12.  = 


(l  +  CT, 

2^ToF;. 

1  +  O-o'fi/ 

inr 

^oll 

i-l  «0i 

For  n  odd 

10-0.05.,T^^ 

T.i  flo* 

For  n  even 

The  series  in  steps  7  and  9  converge  rapidly,  and  three  to  four  terms  are  sufficient  for 
most  purposes.  Using  the  quadratic  formula,  the  i**’  pair  of  complex  pole  values  can  be 
expressed  as 


2 


’0/ 


I  =  1,2,...,/ 


The  zeros  occur  at 


/  =  1,2, 


,/ 
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Appendix  C 


Header  file,  “filtenh” 


/*****♦*******★***♦★*******♦**♦**★★**★★**★♦♦★★♦**♦★*★****★★**♦**★*★ 

♦Purpose;  This  library  defines  the  structure,  functions  and  constants  to 

*  implement  complex  number  manipulation 

*  Author;  Benoit  Gosselin 

*  Date; 

* 

*★***★*******♦★******★★*♦★♦*★★*★★★**★*★•*■**★*★★★★★***★★***★★★★★**♦★/ 

#ifndef  FILTER_H 
tdefine  FILTER_H 

#include  <iostream.h> 
iinclude  <conio.h> 

# include  <math.h> 

Iinclude  <stdlib.h> 

Iinclude  <stdio.h> 

Idefine  PI  3.14159265358979 
Idefine  EPSILON  l.Oe-06 

/*  This  type  implements  a  complex  number  structure*/ 

struct  cnum 

{ 

float  Q; 
float  I; 

}; 

/★*★★★★*♦★**★*******★★★*★★★★★*****★★***★★**★★********♦★***★★★★*★♦** 

★ 

*  Purpose:  This  function  computes  the  real  part  and  the  imaginary 

*  part  of  a  complex  number  expressed  by  r*exp(pi*teta) 

♦  Author:  Benoit  Gosselin 

♦  Date: 

★ 

★  ★***♦****★**★****★*★*****★*★★**■★★****♦***’**★★★*****★***♦♦*■★★**★'*■*/ 


cnum  cfun  (float  c, float  r) 


cnum  nb; 

nb.Q  =  r*cos ( (float) PI*c) ; 
nb.I  =  r*sin( (float) PI*c) ; 
return  nb; 


^★★★★♦♦★♦★★★★★★★★★★★★★★★*******************-**************^********* 

♦ 

*  Purpose:  Generalized  form  for  the  previous  function  to  allow 

*  complex  input 

*  Author:  Benoit  Gosselin 

*  Date: 

* 

*********★★*★****★★****★****★★★*★**★**★*★*****♦★★★★**★★*★★********/ 


cnum  cexp(cnum  c) 

{ 

cnum  nb; 

nb.Q  «  exp (c.Q) *cos (c.I) ; 
nb.I  =  exp(c.Q)*sin(c.l) ; 
return  nb; 

} 

/**★*★***********★★★★★******★★*★★★*★****★*★★********★********■****** 

* 

*  Purpose:  Output  function.  This  function  prints  the  coefficients 

*  computed  for  an  IIR  filter  in  a  text  file. 


Format  of  the  text  file: 


3.000000  //order  of  the  filter 

0.040142  //digital  static  gain 

//numerator *s  coeff  //denominator's  coeff 
0.040142  0.000000  1.000000  0.000000 
0.120425  0.000000  -1.057236  0.000000 
0.120425  0.000000  1.087358  0.000000 
0.040142  0.000000  -0.708990  0.000000 

'1 

first  row  is  for  the  real  part  and  the  second  is  for  imaginary  part 
Benoit  Gosselin 

*  parameters: 

*  cN  :  coefficients  of  the  numerator 

*  cD  :  coefficients  of  the  denominator 

*  N  :  number  of  coefficients  at  the  numerator 

*  M  :  number  of  coefficients  at  the  denominator 

*  order  :  order  of  the  filter 

*  Kd  :  digital  static  gain  of  the  filter 

*  filename  :  name  of  the  file  where  should  be  print  the  coefficients 
**★*★*★♦**★★★★♦*★★**★★**♦*★***★♦♦★*★♦***♦♦★*♦★*★★♦★★★★★*★★★*★♦★**★/ 

void  print_coef f (cnum  cN[],cnum  cD[),int  N,int  M, int  order, float  Kd,char  filename []) 

{ 

int  i; 

FILE  *outputfile, *matoutputfile; 
outputfile  f open  (filename,  "w"); 

//print  an  output  file  in  a  different  format  for  Matlab  uses 
matoutputfile  -  f open ("matcoeffout.txt", "w") ; 
fprintf  (outputfile,  "%f\n**,  (float) order) ; 
fprintf (outputfile, "%f\n", Kd) ; 
cout«endl«"numerator *  s  coefficients " ; 

for(  i  *=  0/i  <  N;  i++) 

{ 

cout«endl«cN[i)  .Q«''  +  i"«cN[i]  .1; 
fprintf (matoutputfile, "%f\n",cN[i] .Q) ; 

fprintf (outputfile, "%f  %f  %f  %f \n", cN [i] .Q, cN [i] . I, cD(i] .Q, cD(i] . I) ; 

} 

for(  i  *=  0/i  <  N;  i++) 

fprintf (matoutputfile, ”%f\n”, cNti) . I) ; 

cout«endl«"denominator *s  coefficients"; 
for(  i  *  0;i  <  M;  i++) 

{ 

cout«endl«cD[i)  .Q«"  +  i"«cD[i]  .1; 
fprintf (matoutputfile, "%f\n",cD[i] .Q) ; 

} 

for(  i  »  0/i  <  M;  i++) 

fprintf (matoutputfile, "%f\n",  cD[i] . I) / 

f close (outputfile) ; 
f close (matoutputfile) / 

} 


the 

Author : 
Date: 


C2 


« 


a 


*  Purpose:  To  multiply  2  polynomials  together 

*  Author:  Benoit  Gosselin 

* 

*  N  and  M  are  the  order  of  polynomials  a  and  b 

*  Date: 

* 


float  *  polym( float  a[],int  N, float  b[l,int  M) 


int  i,j; 

int  length; 

length  =  (N  +  M  +  1); 

i*rj»0; 


float  *poly  «  (float  *)malloc ( (length)  *  sizeof (float) ) ; 

for(i=0;  i<length;i++) 
poly[i]  »  0.0; 

for(i=0;i<N+l;i++) 

{ 

for ( j«0; j<M+l; j++) 

poly[length-i-j-l] +*a[i) *b( j) ; 

) 

return  poly; 


/***★*♦*★*****♦★★*★*★★*★★*★★*★★**★**★★*★★★★*★♦★★****★*♦***★***★★*★* 

* 

*  Purpose:  If  the  real  part  or  the  imaginary  part  of  the  complex 

*  number  is  too  small,  set  it  to  0 

*  Author:  Benoit  Gosselin 

*  Date: 

* 

★*★*♦♦★*★*********♦*******★*★***♦***★★**★***★***★★*♦********★★★*★★/ 


cnum  set_to_zero(cnum  a) 

{ 

if ( f abs (a . Q) <EPSILON) 
a.Q  =  0.0; 

if (fabs(a.I)<EPSILON) 
a. I  *  0.0; 
return  a; 

} 


f 


A 
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*  Purpose:  This  function  computes  the  square  root  of  a  complex  number 

★ 

*  Author:  Benoit  Gosselin 

*  Date: 

* 

★♦★★★*★★*******★★★★**★★********♦★★★*★★*★♦★**♦★***★★★★★****♦♦★♦*★★♦/ 


cnum  sqrtc(cnum  a) 

( 

float  r,  theta/ 

r  »  sqrt(a.Q*a,Q  +  a.I*a.I)/ 

theta  «  atan{a.I/a.Q) ; 

theta  «  theta  /  2.0; 

a.Q  «  sqrt (r) *cos (theta) ; 
a. I  «  sqrt (r) *sin (theta) / 

return  a; 


) 


*  Purpose;  This  function  computes  the  norm  of  a  complex  number 

*  Author:  Benoit  Gosselin 

*  Date: 

* 

**★♦*★★*♦★**★★★★♦★★★★***★★*★*★**♦♦♦♦***♦★*★★*******♦♦**********♦**/ 


float  norm2c  (float  nbQ, float  nbl) 


return  (nbQ*nbQ  +  nbl*nbl); 


/★★★★★★*★★★★*★★***★**★★★★★★*★★♦★★★**★♦♦♦★★******★**♦*************** 

* 

*  Purpose:  This  function  implements  the 

*  multiplication  of  two  complex  numbers 

*  Author;  Benoit  Gosselin 

*  Date: 

* 

*★★★**♦♦*★♦★♦★★★**★★★*♦★*****♦♦**★★★******♦**♦♦***********♦*****♦*/ 


cnum  multc( float  a, float  b, float  c, float  d) 


cnum  temp; 

temp.Q  *  (a  *  c)  -  (b  ♦  d) / 
temp. I  *  (a  *  d)  +  (b  *  c); 
return  temp; 
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/4r  *******************  ****^:<rVlr******#*-********4r**4r*«nk'**4r*<*«-**«*-****** 

* 

*  Purpose:  This  function  implements  the  square  root  of  a  complex  number 

* 

*  Author:  Benoit  Gosselin 

*  Date: 

*  ’ 

cnum  powc  (cnum  a, float  n) 

{ 

float  r,  theta; 
if(n**0) 

{ 

a.Q  «  1; 
a. I  -  0; 
return  a; 

} 

r  -  sqrt(a.Q*a.Q  +  a.I*a.I); 

theta  *  atan(a,I/a.Q) ; 

theta  »  theta  *  n; 

a.Q  *=  pow(r,n)  *cos  (theta) ; 
a. I  *  pow{r,n) *sin (theta) / 

a»set_to_zero (a) ; 
return  a; 


/***★*★**♦***★*****♦***★★**★****♦***★♦**★****♦★*★*********★★******* 

* 

*  Purpose:  This  function  implements  the  division  of  two  complex  numbers 

*  Author:  Benoit  Gosselin 

*  Date: 

* 


cnum  dive (float  NQ, float  NI, float  DQ, float  DI) 


cnum  temp  «  multc (NQ,NI, DQ, -1 .0*DI) ;  //multiplying  by  the  conjugate 

float  norm  *  norm2c (DQ,DI) ; 
if (norm  «=  0) 

{ 

cout«endl«" divide  by  0  in  dive  (norm  *  0)"; 
temp. 0*0/ 
temp. 1*0; 
return  temp; 

} 

else 

{ 

ten^ . Q  *  temp . Q  /  (  norm) ; 
temp. I  *  temp. I  /  (  norm) / 
return  temp; 

) 
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/★*★***★★★★*****★*★★★★*★★**♦★*★★♦*♦*★*♦★★**★*★*★*★★♦★★♦♦**★★★★★*★** 

★ 

*  Purpose:  This  function  in^lements  an  Expansion  recursion  formula  to 

*  obtain  the  coefficients  of  a  polynomial  from  its  roots 

*  Author;  Benoit  Gosselin 

*  Date; 

*  Prameters: 

*  c  :  to  store  the  coefficients  of  the  polynomial 

*  e  :  roots  of  the  polynomial 

*  m  :  order  of  the  polynomial 

*  n  ;  This  indice  decrease  at  each  recursion  loop  from  n-m 

* 


int  coeff  (struct  cnum  ♦c, struct  cnum  etlrint  n,int  m) 

{ 

cnum  num; 

num  -  multc(c[n] .Q,c[n] .I,e(m] .Q,e[ml  .1)  ; 
c(n+l] .Q  *  c[n+l].Q  -  num.Q; 
c(n+l].I  -  c{n+ll.I  ~  num.I; 
if(n  >  0) 

return  coeff (&c[0] , e, n-l,m) / 


else 

} 

#endif 


return  0; 
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